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Simmary 


This  document  is  the  final  report,  for  Phase  TT  of  an  investigation  into 
the  feasibility  and  design  of  a  smart  control  generator  for  the  purpose  of 
automated  image-to-map  and  image-to-image  registration.  This  Phase  II  effort 
is  concerned  with  the  problem  of  extracting  control  information  in  imagery  and 
matching  it  with  corresponding  control  information  in  spatial  databases  in 
order  to  calibrate  the  imagery.  This  problem  is  generically  called  "image-to- 
map  matching".  The  use  of  other  collateral  navigational  data  sources,  such  as 
inertial  navigation  (INS)  and  Global  Positioning  System  (GPS),  is  also 
addressed. 

The  work  describes  the  scenarios  that  will  motivate  the  selection  of  images, 
scenes,  and  collateral  data  and  the  processing  of  this  data  for  generating 
control  information. 
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1.0  Introduction  and  Definition  of  Scenarios 


This  Phase  II  effort  is  concerned  with  the  problem  of  extracting  control 
information  in  imagery  and  matching  it  with  corresponding  control  information 
in  spatial  databases  in  order  to  calibrate  the  imagery.  This  problem  is 
generically  called  "image  to  map  matching."  The  use  of  other  collateral 
navigational  data  sources,  such  as  inertial  navigation  systems  (INS)  euid  Global 
Positioning  System,  is  also  addressed. 

Section  1  describes  the  scenarios  that  will  motivate  the  selection  of  images, 
scenes,  and  collateral  data  and  the  processing  of  this  data  for  generating 
control  information. 

Issues  involving  INS  and  GPS  are  discussed  in  section  2,  along  with  some 
results  of  using  actual  INS  and  GPS  data  for  image  calibration. 

Techniques  for  using  features  in  synthetic  images  generated  from  terrain  data 
to  provide  control  information  is  described  in  section  3. 

Algorithms  for  choosing  and  extracting  spatially  extended  features  in  images 
and  matching  their  descriptions  in  DFAD  data  are  explored  euid  demonstrated  for 
image  control  in  section  4. 

A  discussion  of  the  use  of  image  segments  for  control  information  appears  in 
section  5. 

Techniques  for  calibration  using  satellite  ephemeris  data  are  discussed  in 
section  6. 

A  discussion  of  rule-based  considerations  for  control  information  appears  in 
section  7. 

Spatial  date±iase  issues  are  considered  in  section  8. 

1.1  Definition  of  Mission-Relevant  Scenarios 

Two  tactical  mission  scenarios  drive  the  research  and  test  cases  for  this  Phase 
II  effort,  and  are  distinguished  by  the  presence  or  absence  of  GPS  data.  Both 
scenarios  involve  optical  or  Synthetic  Aperture  Radar  (SAR)  imagery  sensed  from 
aircraft  platforms,  as  distinct  from  satellites.  This  imagery  must  be  properly 
calibrated  to  allow  the  extraction  of  euiy  positional  estimates  of  tactically 
releveint  targets.  ETL's  overall  concept  of  the  various  functions  that  are 
required  for  extracting  tactically  relevant  information  from  imagery  is  shown 
in  Fig.  1.1. 

The  scenes  which  were  selected  for  demonstrating  and  testing  calibration 
concepts  involved  geometrically  extended  features,  variations  in  shading  caused 
by  local  terrain  slopes,  or  shadows  generated  by  occluding  terrain.  The  usual 
ground-based  navigational  aids  such  as  beacons  or  air  traffic  control  radar 
were  assumed  not  to  be  available  for  tactical  mission  scenarios. 
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The  first  tactical  scenario  concerns  the  calibration  of  imagery  using  inertial 
navigation  system  (INS)  data  along  with  a  prior  flight  plaun  euid/or  flight  log. 
The  flight  plan  represents  the  desired  approximate  path  of  the  aircraft.  The 
pilot  will  use  the  output  of  the  INS  2md  visual  sighting  of  landmarks  in 
attenpting  to  execute  this  flight  plain. 

The  INS  positional  data  is  typically  accurate  over  short  time  spans,  but 
suffers  from  drifts  over  longer  periods  due  to  cumulative  integration  of 
errors.  The  heading  information  will  be  assumed  to  be  accurate  to  within  a 
degree,  because  it  is  available  from  other  sources  such  as  con^ss.  Therefore, 
landmark  sightings  can  be  iirportant  for  making  positional  corrections  during 
longer  missions,  and  the  resulting  maneuvers  can  be  recorded  in  a  pilot's  log. 

Imaged  landmarks,  either  point-like  or  extended  features,  and  terrain-induced 
shadows  can  provide  more  accurate  calibration  information.  The  INS  data  and  the 
flight  plan/log  data  provide  approximate  locations  for  the  platform  vhich  can 
be  used  to  initialize  more  precise  position  estimates  using  the  pixel  locations 
of  such  identifieible  features. 

The  other  tactical  scenario  involves  the  use  of  data  from  the  Global 
Positioning  System  (GPS),  along  with  INS  data  for  calibrating  imagery.  This  GPS 
data  caul  be  used  with  the  INS  data  in  a  number  of  ways  in  order  to  improve  the 
platform's  flight  path  estimates.  Because  of  the  high  accuracy  of  GPS,  no 
additional  positional  corrections  using  imagery  should  be  necessary.  Instead, 
the  main  issue  is  how  to  combine  INS  auid  GPS. 

Issues  involving  INS  and  GPS  will  be  discussed  in  section  2.3.1,  and  the  theory 
and  results  of  using  actual  INS  and  GPS  data  for  image  calibration  will  be 
presented  in  sections  2.3.2  and  2.2.3  respectively.  Flight  plan  or  log  data, 
when  available,  can  be  used  to  give  initial  estimates  of  the  flight  path. 

1.2  The  Problem  of  Image-Map  Matching 

In  the  following  discussions,  we  shall  use  the  term  "maps"  in  the  more  general 
sense  of  meaning  map  information  that  is  also  available  in  feature  and  terrain 
date±)ases  as  well  as  in  ordinary  maps  themselves. 

Image-map  matching  is  the  process  of  assigning  a  map  coordinate  to  each  pixel 
in  the  sensed  image.  This  task  can  be  achieved  in  a  number  of  ways.  One 
approach  involves  computing  aui  estimate  of  the  sensor  platform  position  and 
possibly  the  sensor  attitude.  This  problem  is  called  "resection  in  space"  in 
photogrammetric  usage.  Some  other  approaches  are  briefly  mentioned  in  section 
3. 

In  the  past,  resection  in  space  has  been  achieved  with  computational  procedures 
which  relate  individual  pixel  locations  to  ground  control  points  with  known 
geographic  coordinates.  These  procedures  are  based  on  exploiting  point-point 
correspondences.  If  enough  individual  control  points  are  imaged  and  recognized, 
then  such  procedures  are  sufficient.  However,  often  an  imaged  scene  does  not 
contain  the  signatures  of  individually  recognizable  control  points,  but  does 
contain  contours  of  known  objects. 

Little  work  has  been  done  to  generalize  such  resection  in  space  procedures  to 
exploit  contour-contour  correspondences.  This  subject  is  discussed  in  section 

3. 
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The  resection  in  space  problem  requires  the  previous  determination  of  point- 
point  or  contour-contour  correspondences.  This  correspondence  problem  is 
difficult  because  of  a  nxmiber  of  reasons  involving  differences  in  data  types, 
levels  of  detail,  accuracy,  and  data  validity. 

Concerning  the  latter  issue,  it  is  clear  that  maps  are  much  more  static  than 
actual  scenes  in  the  real,  dynamically  changing  world.  For  exanple,  rivers 
meeuider  auid  sometimes  even  dry  up  from  one  season  to  the  next.  Forest  regions 
can  be  burned  or  cut,  urban  areas  can  be  built  up,  and  roads  can  be  built.  Even 
if  scene  content  were  not  unchanged,  as  in  the  aixsve  cases,  images  of  the  same 
scene  can  change  because  of  environmental  effects  such  as  rain  or  snow. 
Clearly,  maps  are  not  updated  at  a  rate  sufficient  to  accomodate  such  changes. 

Maps  and  images  differ  in  the  content  of  their  respective  data.  Maps  generally 
contain  information  of  a  geometric  nature,  whereas  images  are  formed  from 
radiometric  as  well  as  geometric  effects. 

Maps  and  images  also  have  incompatible  data  formats.  Images  are  presented  in 
raster  format  auid  maps  in  vector  format.  Plauiimetric  maps  also  contain  a  great 
deal  of  non-visual  data,  such  as  textual  and  pictorial  annotations,  border 
lines,  etc. 

Real  imaged  scenes  generally  exhibit  wide  varieties  of  rich  detail  at  many 
scales  of  length.  Maps,  on  the  other  hand,  are  abstractions  of  much  less  detail 
at  many  fewer  scales  of  length,  and  contain  information  on  a  more  macroscopic 
level . 

Another  difference  between  the  two  data  types  involves  accuracy  and 
consistency.  Images  can  support  high  measurement  precision  with  high  pixel 
resolution  but  require  collateral  ground  truth  data  to  support  high  accuracy. 
Maps  ceun  be  highly  accurate  in  principle,  but  often  contain  inconsistencies  and 
approximations  which  cam  cause  problems  when  trying  to  merge  these  two  data 
sources . 

For  example,  such  inconsistencies  can  arise  when  performing  updates  to  maps 
using  real  imagery.  A  typical  problem  that  arises,  called  the  "juxtaposition 
problem"  (Goodenough,87] ,  involves  the  outline  of  a  newly  added  feature, 
obtained  from  imagery,  incorrectly  overlaying  the  outline  of  a  feature  already 
on  the  map. 

Such  inconsistencies  would  be  more  tractable  if  there  were  more  descriptive 
information  available  on  the  construction  of  the  map  details.  However,  a  map  is 
an  iconic  abstraction  of  a  real  geographic  region,  and  there  are  in  general  a 
number  of  undocumented  interpolations  and  generalizations  between  the  final 
iconic  representation  and  the  original  data  sources  [McKeown,87 ] . 

Because  of  its  importance  to  numerous  applications,  the  problem  of  image-map 
matching  has  generated  considereible  attention.  In  the  following  section  there 
is  a  discussion  of  some  of  the  work  in  this  area. 
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1.3  Past  and  Recent  Work  in  Image  To  Map  Correspondence 


Because  maps  and  images  have  differing  data  formats  and  data  contents,  any 
atten^ts  at  comparing  and  matching  the  information  between  the  two  types  will 
require  a  common  representation  auid  common  data  content.  Therefore,  the 
structurally  richer,  sensed  image  must  be  processed  so  that  sJsstractions 
compatible  with  map  or  spatial  database  information  are  formed  from  its 
individual  pixel  patterns. 

These  abstractions,  i.e.  features,  from  the  sensed  images  are  then  matched  with 
the  data  in  maps  or  spatial  databases.  The  relevant  issues  then  concern  the 
type  of  features  to  be  used  for  matching,  their  extraction,  their  organization, 
and  the  processes  for  matching  using  these  features.  Some  of  the  particulars  of 
previous  methods  will  be  examined  in  the  course  of  discussing  these  issues. 

In  the  past,  work  on  this  problem  and  on  the  more  general  problem  of 
interpreting  aerial  imagery  was  entirely  concerned  with  the  development  of 
lower-level  image  processing  procedures  for  extraction  auid  matching.  Reviews  of 
such  procedure-oriented  work  can  be  found  in  [Leberl,82],  [Fischler,81] ,  and 
[Medioni,84] .  The  more  general  area  of  model-based  image  analysis  is  thoroughly 
reviewed  in  [Binford,82] ,  amd  more  recent  work  is  selectively  discussed  in 
[Kalvin,86].  Aspects  of  some  of  these  works  are  discussed  in  the  next  sections. 

However,  some  recent  work  has  stemmed  from  the  point  of  view  that  progress  on 
more  general  image  understainding  and  interpretation  problems  cannot  be  achieved 
by  only  working  with  lower  level  data  and  procedures.  Some  sort  of  "higher- 
level  knowledge"  is  required  for  interpreting  images  and  relating  their  content 
to  other  data  sources  such  as  maps. 

Lower-level  processes  based  entirely  on  pixel  luminance  values  can  often  lead 
to  ambiguities  or  errors.  Therefore,  there  may  be  a  requirement  for  a  more 
robust  synthesis  allowing  errors,  inconsistencies,  and  partial  evidence  from 
the  results  of  these  processes.  A  good  example  of  such  a  synthesis  on  a 
procedural  level  can  be  found  in  [Fischler,81 ] ,  and  is  discussed  in  section 
1.3.2. 

Such  a  synthesis  often  depends  strongly  on  the  context  of  the  objects  being 
viewed.  This  context  is  the  "knowledge"  about  real-world  objects  and  their 
relationships.  Recent  work  in  AI  has  focused  on  methods  of  representing  such 
knowledge  and  the  use  of  non-procedural  programming  languages,  such  as  PROLOG, 
for  making  reasoned  conclusions  in  narrowly  focused  domains. 

The  more  general  area  of  aerial  image  "understanding"  is  beyond  the  scope  of 
this  effort.  However,  section  1.3.5  contains  a  discussion  of  some  of  the 
relevant  work  in  knowledge-based  aerial  image  understanding  in  order  to  survey 
some  of  the  ideas  which  can  be  relevant  to  image-map  correspondence. 

1.3.1  Relevant  Features 

Because  of  the  geometric  nature  of  information  stored  in  maps,  the  features  to 
be  abstracted  from  sensed  images  have  been  largely  of  a  geometric  nature  also. 
Intuitively  these  features  involve  extracted  boundaries  of  major  regions  or 
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large  curvilinear  objects,  as  well  as  regions  themselves.  This  intuition  is 
reflected  in  the  practice  of  representing  the  data  stored  in  maps  and  feature 
databases  by  templates  of  Matures,  i.e.  iconic  data.  However,  terrain 
databases  differ  from  these  types  of  data  aund  contain  data  suitaUjle  for  surface 
description. 

The  visible  portions  of  opaque  three  dimensional  objects  correspond  to  points, 
contours,  or  surfaces.  The  two  dimensional  imaged  signatures  of  these  objects 
are  points,  contours,  and  regions.  In  the  past,  the  type  of  feature 
correspondence  that  was  exploited  for  image-map  matching  was  point-point. 
Because  sufficiently  nany  individually  recognizable  control  points  are  not 
always  present  in  imagery,  the  generalization  to  contour-contour  correspondence 
can  become  important.  Contour-contour  mapping  using  DEAD  data  is  discussed  in 
section  4. 

Other  image  signatures  with  definitive  contours  are  terrain-induced  shadows. 
The  use  of  this  type  of  feature  for  image-map  matching  is  described  in  section 

3.2. 


The  extraction  of  contours  depends  on  the  existence  in  the  imagery  of 
relatively  sharp  changes  in  gray-values  rather  than  smooth  chamges.  This 
subject  is  discussed  in  more  detail  in  section  1.3.2. 

On  the  other  hand,  the  use  of  smoothly  changing  image  intensities  does  not  seem 
to  have  been  greatly  exploited  in  image-map  correspondence,  and  has  been 
associated  more  with  robotics  and  industrial  inspection  under  the  heading  of 
"photometric  stereo". 

Photometric  stereo  [Horn, 86]  is  an  example  of  the  usage  of  shading  information 
for  identification  of  industrial  parts.  This  technique  exploits  the 
relationship  between  local  surface  orientations  and  variations  in  image 
intensities. 

Photometric  stereo  is  part  of  a  more  general  subject  area  loosely  called  shape- 
from-shading.  The  use  of  terrain-induced  shading  would  represent  an  option  for 
image-map  matching  if  there  were  a  lack  of  identifiable  points  or  contours 
within  a  sensed  image.  The  use  of  such  terrain-  induced  shading  has  not  been 
explored  for  image-map  matching  to  the  same  extent  as  have  boundary  and  point 
matching. 

However,  one  very  notable  exception  is  described  in  [Horn, 78],  where  optical 
synthetic  images  are  created  using  a  terrain  model.  The  extension  of  this 
method  to  SAR  imagery  is  described  in  section  3. 

1.3.2  Extraction  of  Relevant  Features 

Clearly,  the  extraction  of  features  in  the  sensed  imagery  should  not  proceed 
independently  of  the  knowledge  represented  in  the  maps  and  databases,  and  any 
approximate  knowledge  of  sensor  platform  position.  Unguided  segmentation  based 
solely  on  image  luminance  values  has  long  been  known  to  lead  to  problems 
because  object  and  intensity  boundaries  are  not  always  the  same.  The  use  of 
this  collateral  map  knowledge  is  loosely  referred  to  as  "map-guided 
segmentation",  a  term  coined  in  [ Barrows, 77 ] . 
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The  approaches  generally  taken  involve  either  boundary  fragment  extraction 
followed  by  some  form  of  boundary  matching  or  template  matching. 

Binary  masks  were  used  in  [Kropatsch,81]  for  feature  extraction  using  templates 
of  objects  known  to  have  been  in  the  vicinity  of  the  imaged  area.  These  masks 
were  apparently  not  modified  to  take  into  account  ciny  projective  effects, 
however. 

The  MAPS  system  [McKeown,87]  also  uses  a  priori  location  information  to  infer 
likely  features  present  in  the  sensed  image.  It  then  projects  models  of  these 
features  into  the  image  using  the  camera  model,  to  aid  in  the  extraction  of 
these  features  eind  to  make  displacement  calculations. 

Two  categories  of  edge  extraction  operators  were  used  in  [ Fischler ,81 ]  for 
extracting  roads  in  optical  images.  Type  1  operators  were  classified  as  those 
which  never  incorrectly  classify,  but  which  sometimes  miss  correct  instances. 
Type  2  operators  were  those  which  sometimes  misclassified,  but  which  accurately 
measure  true  insteinces. 

The  output  of  these  operator  types  is  combined  using  a  cost  function  to  give  a 
cost  measure  of  a  road  going  through  every  pixel  in  a  some  area.  This  cost 
faction  can  incorporate  some  prior  knowledge  about  road  curvature  and  other 
parameters  of  road  structure  as  weighting  terms  in  this  cost  function.  The 
extracted  road  is  obtained  as  a  least  cost  path  through  the  connected  graph  of 
possible  road  locations.  This  least  cost  path  is  calculated  using  F*,  an 
efficient  modification  of  the  heuristic  graph  search  algorithm  A*  [Pearl, 84]. 

The  importance  of  the  method  used  in  [Fischler, 81]  is  the  use  of  multiple, 
possibly  contradicting,  local  sources  of  evidence,  and  using  them  to  estimate  a 
more  global  feature  based  on  a  minimization  criterion. 

1.3.3  Feature  Matching  Techniques 

We  will  classify  feature  matching  techniques  as  region-based  or  boundary-based, 
and  these  are  discussed  in  sections  1.3. 3.1  and  1.3. 3. 2  respectively.  In  either 
case,  matching  requires  some  analytic  metric  for  making  ranked  comparisons  of 
good,  poor,  and  ambiguous  matches.  The  challenge  is  to  find  metrics  which  are 
robust  with  respect  to  a  wide  variety  of  illumination  and  viewing  conditions. 

For  the  image-map  matching  problem,  the  approach  is  to  infer  some  relationship 
of  extracted  features  within  the  image  and  find  the  corresponding  relationship 
within  the  map. 

1.3. 3.1  Region-Based  Matching 

Template  matching  using  any  of  the  various  forms  of  the  correlation  metric  has 
historically  been  a  useful  method  for  image-image  matching.  This  method  of 
matching  is  generally  region-based  rather  than  boundary- based,  although  it  can 
also  be  used  with  some  degree  of  success  for  matching  high-pass  filtered  data. 
Another  area-based  metric  is  the  sum  of  absolute  differences  [Barnea,72). 

One  classification  scheme  of  image  translation  estimation  methods  appears  in 
[Huang, 81 1,  and  consists  of:  Fourier-based,  matching,  and  differential. 
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Differential  methods  apply  to  image  pairs  with  small  amoionts  of  relative 
translation,  such  as  successive  video  frames.  Such  methods  are  based  on 
Taylor's  expeinsion  for  two  vari2d3les  truncated  to  linear  terms.  Obviously,  such 
cin  approach  is  not  applicable  to  this  registration  problem,  and  will  not  be 
discussed  here. 

Matching-type  methods  refer  to  the  use  of  correlation-based  methods  as  a 
measure  of  similarity  for  trial  areas  of  overlap  between  image  pairs.  The  use 
of  such  methods  will  be  discussed  in  sections  3.2,  3.3,  4.  and  5. 

Fcurier-based  methods  refer  to  an  explicit  use  of  the  Fourier  Shift  Theorem 
[Champeney,73 ] ,  which  states  that  if  two  integrable  continuous  fionctions  f(x,y) 
and  g(x,y)  are  related  as:  g(x,y)  -  f(x  +  /ix,  y  +  Ay)  then  their  Fourier 
Trcinsforms  are  related  as: 

F(u,v)  -  G(u,v)  exp{  -j  2II  (u  Ax  +  V  Ay)) 

In  practice,  such  methods  reduce  to  phase  correlation. 

A  classical  method  for  registering  images  is  to  successively  register 
corresponding  patches  within  the  two  images  using  an  analytic  metric  for 
comparison.  For  each  such  patch  in  one  image,  candidate  trial  patches  in  the 
second  image  are  compared  using  this  metric.  The  trial  patch  in  the  second 
image  which  optimizes  the  metric  is  chosen  as  the  "matching"  patch  for  the 
given  patch  in  the  first  image. 

Historically,  the  normalized  correlation  coefficient  has  been  the  preferred 
metric,  although  other  less  robust  versions  of  correlation-type  metrics  have 
been  used  because  of  computational  advantages,  such  as  non-no rmali zed 
correlation,  or  variations  on  the  sura  of  absolute  difference  metrics.  These 
three  are  given  below: 

Normalized  Correlation: 

E  Z  F  (j,k)  F  (j-m,k-n) 

j  -  1  k  -  1 

R(m,n)  -  - 

[Z  Z  F  { j,k)  ]  [ Z  Z  F  ( j-m,k-n) ] 

3-lk-l 

Correlation: 

R(m,n)  -  Z  Z  F  {j,k)  F  (j-m,k-n) 

:-i  k-i 

Sum  of  Absolute  Differences: 


R(m,n)  =  Z  Z  IF  (j,k)  -  F  (j-m,k-n)‘ 


Normalized  correlation  is  the  preferred  method  vdienever  computational  cost  is 
not  an  issue.  The  value  of  the  normalized  correlation  coefficient  is  bounded 
absolutely,  and  always  lies  between  -1  and  +1.  Therefore,  2dDSolute  values  of 
this  metric  close  to  +1  achieved  for  trial  offsets  are  close  to  local  naxima  in 
the  correlation  surface.  The  normalized  metric  cein  allow  for  a  constant 
multiplier  between  the  gray-values  for  corresponding  pixels  in  the  two  images. 
Also,  the  normalized  correlation  coefficient  theoretically  achieves  its  maximum 
value,  in  the  absense  of  noise  auid  nonlinearities,  at  the  correct  offset 
between  the  two  images. 

This  is  not  always  the  case  for  non-normal i zed  correlation  algoritiims,  which 
can  achieve  high  correlation  values  at  some  incorrect  offsets  simply  because  of 
high  gray-values  at  certain  localities.  The  sum  of  absolute  values  of 
differences  generally  performs  better,  in  the  sense  of  acceptable  accuracy, 
than  Lhe  non-normal i zed  correlation  approach,  but  also  not  as  well  as 
normalized  correlation  [ Svedlow, 78 ] . 

Correlation  is  computationally  expensive.  One  route  toward  reducing  the  amount 
of  computation  has  been  to  perform  the  equivalent  operation  in  the  Fourier 
domain,  using  the  Fourier  Convolution  Theorem  [Champeney,73 ] .  The  basis  of 
computational  savings  using  the  FFT  is  that  the  latter  requires  a  computation 
asymptotically  proportional  to  N  log  N  rather  than  N^ ,  where  N  is  the  dimension 
of  a  square  image.  For  relatively  large  images,  there  cein  be  consideredDle 
savings  using  the  FFT  approach. 

However,  this  approach  does  not  deal  with  the  observation  that  the 
determination  whether  a  certain  offset  is  very  incorrect  should  somehow  require 
less  computation  than  the  determination  of  a  correct  offset.  This  view  is  at 
the  root  of  the  idea  of  recursively  searching  for  the  correct  offset  in  a 
"pyramid"  of  increasing-resolution  versions  of  the  image  pairs  [Rosenfeld,84  ] . 

This  point-of-view  is  also  the  basis  of  the  sequential  approach  using  the  sum 
of  a±)solute  differences  metric  in  [Barnea,72].  Here,  for  any  trial  offset, 
whenever  a  pre-determined  threshold  has  been  exceeded  before  the  entire  sum  has 
been  evaluated,  the  summation  is  suspended  and  a  new  trial  offset  is  evaluated. 

Hybrid  algorithms  using  this  approach  for  a  rough  estimate  of  offset,  followed 
by  a  version  of  normalized  correlation  using  statistical  pre-processing  have 
been  suggested  [Pratt, 73].  Clearly,  this  idea  can  be  generalized  to  the  use  of 
other  robust  metrics  following  the  initial  rough  estimate. 

Another  problem  with  correlation  measures  in  general,  including  normalized 
correlation,  includes  the  broad,  flat  nature  of  the  peak  regions  in  the 
correlation  surface.  This  characteristic  has  advantages  and  disadvantages.  An 
advantage  of  broad  peaks  is  their  "pull-in"  range  for  search  techniques  that 
don't  sample  everywhere  or  that  employ  averaging,  as  in  multi-resolution 
processing  [ Rosenfeld,84 ] . 

A  disadvantage  is  that  registration  accuracy  can  be  compromised  whenever  peaks 
arc  not  sharp,  and  smaller  perturbations  due  to  noise  can  potentially  have 
greater  effects  on  accuracy. 
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This  broad  peak  characteristic  of  correlation  techniques  occurs  because  spatial 
relationships  in  images  are  ignored  for  the  most  part.  Instead,  the  values  of 
the  correlation  metric  are  really  related  to  energies  of  the  broad  areas  within 
the  images,  i.e.  the  energy  content  of  the  low-frequency  portions  of  the 
images.  Since  phase  shifts  of  these  lower  frequencies  cam  be  relatively  large 
compared  to  the  pixel  resolution  and  still  be  relatively  small  conpared  to  the 
corresponding  low-frequency  wavelengths,  there  is  a  lack  of  sensitivity  in 
correlation  metrics  to  phase  shifts.  Therefore,  correlation  surfaces  tend  to 
have  broad  peaks. 

One  approach  toward  solving  this  problem  has  involved  the  preferential  use  of 
phase  information  in  the  images.  The  idea  here  is  that  the  Fourier  phase 
content  of  the  images  contains  more  accurate  information  than  the  low- 
frequencies  which  dominate  correlation  metrics.  This  is  the  basis  for  phase 
correlation  methods  [Kuglin,75,79] ,  [Pearson, 77] ,  [DeCastro,87 ] . 

Using  the  previous  example  of  the  functions  f(x,y)  and  g(x,y)  related  by 
translations  Ax  and  Ay,  the  inverse  Fourier  transform  of  {F(u,v)  /  G(u,v)} 
is  just  the  inverse  transform  of  the  phase  term,  and  thus  is  equal  to  the  Dirac 
Delta  Distribution  evaluated  at  Ax  and  Ay,  i.e.  5(Ax,  Ay) 

Now,  because  of  the  effects  of  sampling  and  finite  image  sizes,  sidelobes  occur 
in  addition  to  a  main  peak.  Therefore,  in  practice,  this  technique  reduces  to 
correlation  in  the  Fourier  phase  domain. 

Such  methods  are  also  capable  of  subpixel  accuracy  with  the  use  of 
interpolation,  but  suffer  from  the  problem  of  potentially  high  sidelobes.  Such 
a  correlation  surface,  with  sharp  peaks  and  high  sidelobes,  does  not  lend 
itself  easily  to  hierarchical  processing  with  reduced  resolutions  because  of 
the  narrow  "pull-in"  rcinge  of  the  main  lobe.  However,  used  in  conjunction  with 
other  methods  which  can  acquire  the  main  lobe,  phase  correlation  can  be  a 
useful  method  for  refining  initial  offset  estimates  to  subpixel  accuracy. 

Another  approach  to  rectify  this  broad  peak  problem  has  been  to  concentrate  on 
the  edge  content  in  the  images.  This  approach  has  both  intuitive  appeal  as  well 
as  some  theoretical  justi  f  i  r^t-ion. 

Intuitively,  it  would  seem  that  the  "meaningful"  information  in  an  image  lies 
at  the  locations  of  large,  structurally  significant  contours  and  boundaries. 
Particularly  in  SAR  images,  the  smaller  edges  are  more  often  due  to  noise, 
speckle,  and  imaging  effects.  However,  the  larger  edges  are  due  to  terrain  and 
thematic  value  changes.  It  would  seem  that  it  is  these  edges  and  boundaries 
that  should  be  involved  in  a  registration  scheme. 

This  notion  has  been  explored  on  a  more  analytical  basis  from  tho  stanHpoint  of 
"optimally"  filtering  the  images  as  a  pre-processing  step  prior  to 
registration.  One  approach  has  been  to  decorrelate  the  images  by  applying  a 
"■whitening"  filter  [Pratt, 73],  [Svedlow,78 ] .  These  methods  essentially  differ 
in  their  assumptions  concerning  image  structure  and  statistical  properties.  The 
conclusions  of  these  works  point  toward  the  use  of  pre-processing  filters  which 
can  be  approximated,  'under  certain  assumptions,  by  gradient  filters 
[ Svedlow, '8  I ,  or  Laplacian  filters  [Pratt, 73]. 
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Of  course,  there  are  some  problems  associated  with  this  approach  also.  The 
high-frequency  edge  content  of  an  image  is  relatively  small  compared  to  the 
total  area  of  an  image.  Edges  can  be  broken  up  slightly  differently  in  two 
images  which  otherwise  contain  no  other  perturbations.  Edges  may  also  have 
slight  variations  in  thickness.  Therefore,  such  perturbations  can  lead  to 
misalignment  sensitivities  for  algorithms,  like  edge  correlation,  v^ich  examine 
the  degree  of  "match"  in  overlapping  edge  images.  One  approach  is  to  condition 
the  edges  in  both  images.  This  would  include  normalizing  their  intensities  and 
broadening  them.  However,  broadening  edges  can  lead  to  reduction  of 
registration  accuracy. 

These  considerations  mentioned  above  suggest  that  correlation-based  matching 
involves  consideraible  searching  and  is  computationally  expensive,  unless  the 
hierarchical  "coarse-to-fine"  approach  previously  discussed  is  used.  The 
capabilities  of  correlation  algorithms  for  achieving  high  accuracies  often 
seems  to  require  some  form  of  pre-filtering  to  enhance  edge  content.  However, 
certain  potential  instaibilities  are  involved  with  the  use  of  edge  correlation. 

In  [ Kropatsch, 81 ]  normalized  correlation  with  a  binary  mask  matrix  was  used  for 
both  feature  extraction  as  well  as  feature  matching  in  image-map  matching.  The 
use  of  correlation,  however,  assumes  that  there  is  little  nonlinear  geometric 
or  radiometric  distortion  between  the  template  and  the  actual  imaged  feature 
instance  (Lahart,70]. 

Synthetic  images  corresponding  to  arbitrary  viewing  scenarios  can  be  created 
using  a  terrain  database.  The  available  information  in  these  synthetic  images 
consists  of  shading  corresponding  to  the  interaction  of  the  macro-relief  of  the 
terrain  model  and  the  viewing  geometry. 

Included  in  these  shaded  regions  would  be  terrain-induced  shadows.  The  shadows, 
although  corrupted  by  noise,  are  more  distinctive  and  should  not  involve  other 
types  of  micro-structure.  However,  shadows  are  potentially  more  susceptible  to 
certain  instabilities.  The  use  of  shadows  for  image-map  matching  is  discussed 
in  section  3. 

The  use  of  synthetic  images  created  from  terrain  databases  for  optical  imagery 
is  described  in  [Horn, 78].  Two  similar  matching  techniques  were  used.  The 
first  used  the  normalized  correlation  metric  which  normalizes  the  correlation 
coefficient  by  the  stcindard  deviations  of  the  two  images.  The  second  method 
consisted  of  a  modification  of  the  normalizing  term  using  the  arithmetic  mean 
[Moravec, 77 ] ,  which  was  easier  to  compute.  The  registration  results  were 
insensitive  to  the  choice  of  normalizing  terms. 

As  a  variation  on  this  theme,  the  use  of  synthetic  images  created  from 
planimetric  maps  was  also  attempted  in  [Triendl,81]  but  specifically  rejected 
as  a  general  method  in  [Leberl,82].  The  reason  given  for  rejection  was  that 
maps  contain  a  great  deal  of  non-image  data.  Also,  it  was  clear  that  matching 
with  the  sensed  images  would  also  require  considerable  filtering  to  remove  the 
smaller  levels  of  micro-structure  present  in  imagery  which  are  never  present  in 
maps . 


1.3. 3.2  Boundary-Based  Matching 


One  methcxl  of  boundary  matching  in  binarized  images  which  has  achieved  some 
success  is  "chamfer  matching"  [Barrows,  78].  This  method  uses  a  distauice  array 
between  features  in  patches  being  matched,  and  estimates  the  tramslations 
between  patches  by  searching  for  those  offsets  which  minimize  the  distance 
array  sums.  This  technique  generally  requires  cleaner  extraction  of  edges  than 
does  edge  correlation.  As  is  the  case  with  all  edge-based  matching  techniques, 
this  metric  is  more  sensitive  to  perturbations  perpendicular  to  lines  in  the 
image  than  parallel  to  such  lines. 

One  method  developed  at  VEXCEL  [ McConnell, 87 ]  performs  matching  for  translation 
estimation  by  creating  binarized  images  from  Marr-Hildreth  operator  zero 
crossings.  This  method  has  the  advantage  of  accurate  computation  of  edge 
content  using  this  second-order  operator,  as  well  as  the  stability  of  regional 
matching.  This  method  is  especially  well-suited  for  registering  opposite-side 
satellite  SAR  images. 

Another  method  was  developed  at  VEXCEL  for  registering  ice  floes  in  arctic  SAR 
imagery.  This  method  performed  boundary  matching  using  dynamic  programming,  and 
reduced  a  two-dimensional  search  problem  to  one  involving  one-dimension.  This 
method  employed  the  psi-s  representation  of  closed  contours  [Ballard, 82 ] ,  which 
allowed  a  convenient  representation  for  translation  and  rotation.  The 
particular  implementation  of  the  dynamic  programming  algorithm  [Sankoff,83]  was 
that  used  for  string  search  and  other  sequence  comparisons.  The  assumption  was 
that  the  ice  floes  are  rigidly  rotated  and  translated,  but  some  very  localized 
distortions  such  as  expansion,  contraction,  insertion  and  deletion  could  be 
tolerated. 

Graph  matching  is  a  method  which  incorporates  topological  relationships  into 
the  matching  process  without  undue  emphasis  on  "exact"  metric  correspondence. 
This  method  has  found  considerable  use  in  image-image  matching,  for  exeunple 
[Price,  82],  [Ballard, 82] ,  [Nevatia,  82],  [Shapiro, 81] ,  [Ayache,87].  Such 
problems  involve  subgraph  isomorphism  and  can  potentially  be  complex.  Although 
graph  matching  problems  belong  to  the  worst-case  computationally  intractable 
class  NP-complete  [Aho,74],  the  use  of  heuristics  to  reduce  the  search  has  been 
successful . 

In  fact,  all  computationally  successful  graph  matching  algorithms  must  somehow 
use  some  means  of  distinguishing  salient  features  in  order  to  avoid  the  problem 
of  combinatorial  explosion  when  searching. 

For  example,  in  [Medioni ,84 ] ,  graph  matching  was  used  for  both  image-image  and 
image-map  matching.  The  computational  complexity  of  matching  was  reduced  by  the 
use  of  a  "coarse  to  fine"  matching  strategy  which  extracted  isolated,  long 
edges  first,  and  matched  these  to  the  model  at  low  resolution  using  relaxation. 

Other  boundary-oriented  graph  matching  techniques  have  emerged  from  the  field 
of  robot  vision.  Typical  of  the  problems  encountered  in  this  field  is  to 
correctly  identify  the  individual  parts  in  an  imaged  pile  of  parts  which  may  be 
overlapping  and  are  partially  obscured  by  each  other.  Since  the  surfaces  of 
these  parts  are  often  smooth,  the  boundaries  contain  most  of  the  identifying 
information  (although  see  section  1.3.1  for  a  brief  discussion  of  photometric 
stereo ) . 
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An  example  of  such  work  is  (Holies, 82).  Significant  local  features  are 
identified,  such  as  corners  and  holes.  These  features  are  then  clustered,  and 
matching  proceeds  "cluster  to  cluster".  The  algorithm  attempts  to  find  the  most 
significant  cluster.  A  graph  matching  formulation  is  created  wherein  the  nodes 
are  matches  of  image  and  model  features,  and  edges  are  pair-wise  assignments 
between  nodes.  The  match  problem  is  then  equivalent  to  searching  for  maximal 
cliques,  \diich  has  complexity  NP-complete.  This  method  assumes  that  features 
will  be  clustered  closely  together. 

Another  method  is  that  found  in  [Ayache,84].  Here,  polygonal  approximations  of 
extracted  boundaries  and  model  sides  are  matched  using  a  strategy  of  matching 
the  longest  sides  first,  as  is  done  in  [Medioni ,84 ) .  Lengths  of  sides  and 
corner  angles  are  used  in  measures  of  "compatibility".  This  method  also 
attempts  to  generate  "hypotheses"  and  continues  until  a  sufficient  number  of 
hypotheses  are  evaluated  and  a  good  match  score  has  been  obtained. 

The  method  of  [Turney, 85]  divides  a  template  of  length  n  into  n/2  subtemplates 
of  length  h.  Ever^*  other  pixel  on  the  boundary  of  the  template  starts  a 
subtemplate.  Eacn  subtemplate  of  each  model  object  is  then  compared  to  the 
extracted  ooject  boundary  in  the  sensed  image. 

The  boundaries  are  represented  using  a  variation  of  the  psi-s  representation 
mentioned  previously.  The  matching  metric  is  the  least-squares  measure  of  the 
differences.  The  method  attempts  to  reduce  the  combinatorics  of  the  comparisons 
by  selecting  "most  salient"  subtemplates.  The  successful  matching  of  such 
salient  subtemplates  is  weighted  more  than  other  matches.  Again,  this  algorithm 
is  slow. 

The  spirit  of  this  last  method,  however,  is  carried  out  with  significant 
computational  improvements  in  [Schwarz, 85] .  Moreover,  [Kalvin,86]  continues 
computational  improvements  on  this  approach  by  employing  the  clever  idea  of 
"geometric  hashing".  This  hashing  is  accomplished  by  the  use  of  a  "footprint" 
of  an  object's  boundary.  Only  a  small  number  of  object  footprints  are  retrieved 
for  potential  matching  checks.  This  number  of  retrieved  footprints  does  not 
depend  strongly  on  the  number  of  objects  in  the  model  database. 

A  footprint  is  a  crude  geometric  characterization  of  an  object  boundary  using  5 
dimensions.  A  footprint  is  generated  by  a  mapping  that  is  not  necessarily  1-1, 
but  satisfies  invariance  under  translation  and  rotation,  continuity  is 
essentially  preserved  by  having  locally  similar  objects  into  similar 
footprints.  The  five  dimensional  representation  of  each  point  of  the  object's 
boundary  involves  the  first  four  Fourier  coefficients  of  the  boundary  and  the 
turning  angle  of  a  polygonal  approximation  at  that  boundary  point. 

The  hashing  occurs  because  5-D  space  is  divided  into  hypercubes  of  fixed  si:^e, 
and  for  each  hypercube  there  is  a  list  of  all  models  whose  footprints  lie  in 
that  hypercube.  In  this  way,  searching  for  potential  matches  of  Icoundary 
segments  is  considerably  reduced. 
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1.3.4  Knowledge-Based  Aerial  Image  /Analysis 

One  of  the  tenets  of  this  point  of  view  is  that  a  single  pass  over  the  data 
with  any  particular  extraction  algorithm  will  generally  not  be  conclusive. 
Instead,  wliat  is  required  is  an  iterative  process  that  examines  results  from 
one  or  more  procedures,  formulates  hypotheses,  and  gradually  eliminates 
competing  hypotheses  by  a  "convergence  of  evidence". 

In  particular,  some  importeint  considerations  for  matching  processes  are  the  use 
of  domain-specific  collateral  knowledge,  partial  matching,  hypothesis 
formation,  the  use  and  propagation  of  constraints,  selective  attention  and 
searching  strategies,  accumulation  of  confirming  evidence,  resolution  of 
conflicting  evidence,  and  efficiency. 
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2.0  Generation  of  Control  Information  Using  INS/GPS 

This  section  discusses  the  issues  associated  with  the  use  of  dead-reckoning 
aids  such  as  inertial-  navigation  systems  (INS)  as  well  as  the  Global 
Positioning  System  (GPS),  also  sometimes  called  NAVSTAR,  for  calibrating 
imagery. 

The  issues  associated  with  INS  and  GPS  separately  are  discussed  in  sections  2.1 
and  2.2  respectively.  The  use  of  INS  and  GPS  together  is  discussed  in  section 

2.3.1  along  with  some  results  in  the  open  literature.  An  example  of  the  use  of 
filtered  INS/GPS  data  for  calibrating  actual  detected  SAR  imagery  is  presented 
in  section  2.3.2.  Finally,  a  brief  discussion  on  the  use  of  INS/GPS  data  for 
making  real-time  phase  adjustments  to  raw,  complex  SAR  data  is  given  in  section 
2.3.3. 

2.1  Inertial  Navigation  Systems 

An  inertial  navigation  system  (INS)  contains  an  inertial  measuring  unit  (IMU) 
auid  equipment  for  the  stadjilization  eund  processing  of  sensor  outputs  from  the 
IMU.  The  IMU  provides  sensor  instruments  for  the  three-dimensional  monitoring 
of  absolute  rotation  cind  non-gravitational  translational  acceleration.  The 
latter  is  measured  relative  to  a  stationary  point  with  respect  to  an  inertial 
coordinate  system  or  a  point  moving  with  constant  velocity  with  respect  to  this 
coordinate  system. 

The  instnoments  in  an  IMU  are  accelerometers  for  measuring  linear  accelerations 
and  gyroscopes  for  measuring  rotational  rates.  The  IMU  is  either  implemented  as 
a  gimballed  or  strapdown  system.  A  gimballed  system  is  mounted  on  a  servo- 
driven  platform,  whereas  a  strapdown  system  is  attached  directly  to  the  vehicle 
relative  to  a  vehicle-based  coordinate  system.  Gimballed  systems  are  more 
expensive  than  strapdown,  but  are  more  accurate  over  longer  time  periods  on  the 
scale  of  hours.  For  flight  times  on  the  scale  of  minutes,  strapdown  systems  are 
generally  adequate.  Errors  also  propagate  differently  for  both  types  of 
systems . 

A  gimballed  system  is  siabjected  only  to  relatively  small,  short  angular  rates 
because  the  servos  counteract  these  rotations  in  order  to  preserve  the 
"leveling"  of  the  system.  Moreover,  only  small  torquers  are  required  since  the 
incremental  rotations  are  small. 


On  the  other  hand,  the  strapdown  system  is  subjected  to  all  the  vehicle 
rotations.  The  gyro  rebalancing  signals  must  be  computed,  requiring  larger 
torquers.  The  accuracy  of  integrating  these  rebalancing  signals  is  not  as  great 
as  is  obtained  by  the  incre»",ental  stabilization  of  the  gimballed  platform. 
Likewise,  the  accelerometer  directions  are  incrementally  held  along  the 
reference  axes  in  a  gimballed  system,  whereas  these  accelerometer  outputs  must 
be  recomputed  onto  reference  directions  in  a  strapdown  system. 

2.1.1  Gyroscope  Technologies 

The  most  complex  component  of  an  IMU  is  the  gyroscope.  In  the  past,  the 
mechanical  gyro  has  been  used  in  an  INS.  The  mechanical  gyro  is,  of  course, 
presently  the  best  understood  and  perhaps  still  the  most  reliable  type. 
However,  it  also  has  the  greatest  power  and  weight  requirements  and  involves 
mechanically  moving  parts. 


However,  newer  types  of  gyros  based  on  laser  technology  have  been  in  various 
stages  of  development.  Such  technologies  do  not  depend  on  moving  parts  and 
friction,  and  include  the  ring  laser  gyro,  the  fiber  optic  gyro,  and  the 
hemispherical  resonant  gyro  (Lundberg,87] . 

The  ring  laser  gyro  is  for  strapdown  systems,  and  its  operating  principles  are 
based  on  the  interference  pattern  created  by  two  planar,  oppositely  circulating 
laser  signals.  Mirrors  are  used  to  control  the  paths  of  these  signals.  By 
measuring  phase  shifts  of  the  interference  pattern,  rotation  about  an  axis 
perpendicular  to  this  plane  can  be  estimated.  The  corresponding  power  and 
weight  requirements  are  less  than  for  mechanical  gyros.  However,  the 
observability  of  low  rates  of  rotation  is  compromised  by  phase  lock,  and 
requires  mechanical  dithering  to  counteract  [Lmdberg,87] .  Such  a  system  has 
been  in  general  service  with  some  airline  aircraft  since  1982  [Divakaruni,87] . 

The  fiber  optic  gyro  is  based  on  the  same  physical  principles  as  the  ring  laser 
gyro,  but  does  not  use  mirrors  to  control  the  laser  signal  paths.  Instead,  a 
fiber  optic  coil  is  used.  No  phase  lock  problem  occurs  at  low  rotation  rates 
fLundberg,87] . 

Finally,  the  hemispherical  resonanat  gyro  is  based  on  the  principle  that 
standing,  resonant  waves  on  a  hemisphere  do  not  rotate  at  the  same  rate  as  the 
hemisphere.  Some  phase  locking  occurs  at  lower  rotation  rates,  but  is  corrected 
using  closed  loop  control  [Lundberg,87] ,  [Loper,87]. 

2.1.2  INS  Error  Sources 

The  error  sources  in  an  INS  include  bias  errors  and  other  degradations  in  the 
accelerometers  and  gyros.  The  following  discussion  summarizes  some  arguments 
in  [Farrell,  76],  in  order  to  give  a  very  simplified  overview  of  such  errors. 

Bias  errors  due  to  incorrect  initializations  give  rise  to  oscillations  at  the 
Schuler  frequency.  To  see  this,  suppose  a  plane  is  travelling  on  an  eastbound 
equatorial  path  with  constant  altitude.  The  platform  must  be  kept  level  with  a 
leveling  command  of  Vj./B.  However,  an  initial  nonzero  uncertainty  of: 

Vj.  would  overdrive  the  platform  at  a  rate  of: 

ir,  -  -  V,/R 

where  Vj.  -  Vj.  -  Vj. 

Vj.  is  estimate  of  Vg 

This  incorrect  tilt  would  result  in  aerodynamic  lift,  causing  a  retarding 
acceleration: 

Vg  -  Vg  -  Vg  =  g%  (^'9  given  below) 

Differentiating  and  combining  these  equations  leads  to: 

Vg  +  (g/R)  Vg  =  ij/j,  +  (g/R)  tj,  =  0 
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The  solutions  to  this  linear,  second  order  differential  equation  contains 
sinusoids  at  the  Schuler  frequency  W-(g/R)^/^  rad/sec.  More  generally,  it  can 
be  shovm  that: 

cos  Wt  +  +  (n,j/W))  sin  Wt 

+  Rd  (1-cos  Wt)  +  RD 

COS  Wt  +  ((d/W)  -  (Vj.,^)/u))  sin  Wt 
-(n.2/9)  +  (Wj  >\3/WM)  (1  -  cos  wt) 

where:  Vj, ,  V^.^)  are  velocity  error  and  initial  velocity  error,  respectively 

in  East  direction 

W  is  Schuler  frequency 

g  is  gravitational  acceleration 

R  is  nominal  distance  from  center  of  Earth 

u  -  /g/R 

'^A(o)  respectively,  the  inertial  system  attitude 

uncertainty  vector  with  respect  to  geographic  reference  system,  and 
the  initial  uncertainties  with  respect  to  North,  and  North  x  East 

n, 2  is  the  East  conponent  of  the  acceleration  error 

'  ^2'  *^3  components  of  the  rotational  drift  rate 

Wj ,  Wj ,  Wj  are  the  components  of  the  euigular  drift  rate  of  the 
reference  frame 

d  -  n„i  -  Wj 

D  is  the  time-varying  drift  rate 
-  w^n^j  ( (sinWt/W)-t) 

Now,  except  for  the  last  term  involving  azimuth  drift  rate,  the  velocity  error 
consists  of  sinusoidal  oscillations  at  (W/2n)  Hz,  the  Schuler  frequency. 

Interestingly,  the  tilt  rate  and  accelerometer  biases  do  not  give  rise  to 
continuously  increasing  tilt  or  velocity  errors.  The  velocity  error  produces  a 
positional  error  in  the  East  direction  of: 

K  -  (sin  Wt/W)+R(v|^,^,+n,2/g)(l  cost) 

+RD(t-(sin  Wt/W)  )+R  Wj  n^,3  ( (1-cos  Wt/w^ )  -  ) 
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Now,  for  time  durations  of  <  (2II)/(10W),  the  following  approximations  are 
valid: 

sin  Wt/W  =  t  -  (1/6)  (Wt)2  t 
1-cos  Wt  =  (1/2)  (Wt)2 

so:  Xj  =  Xj,,,+VE,„,+l/2(gtj^,„,+n,2  )t2+Rd(st)2t 

This  last  expression  is  now  decomposed  in  terms  of  eui  initial  position  error,  a 
linear  term  proportional  to  the  initial  velocity  error,  a  quadratic  term 
proportional  to  initial  acceleration  uncertainty,  and  a  drift  rate  of  (1/6) 
Rd(Wt).  This  last  term  is  about  an  order  of  magnitude  less  than  d  for  time 
durations  of  about  1/6  the  Schuler  period,  which  is  about  84.4  min  [Meirovitch, 
70].  This  meeuis  that  a  relatively  high  drift  rate  could  be  tolerated  for  such 
short  time  periods. 

Errors  also  arise  from  uncompensated  gravitational  uncertainties,  which  are 
equivalent  to  time-varying  accelerometer  errors.  Other  errors  also  arise  from 
non-const2uit  degradations  of  the  components,  such  as  the  accelerometers  and 
gyros. 

Narrow  beuid  noise  auid  these  components  can  also  lead  to  bias  errors  through  a 
process  called  rectification  (not  to  be  confused  with  the  photogrammetric 
term),  as  well  as  to  random  rotational  drift  rate  and  acceleration  noise. 
Rectification  is  the  process  of  creating  constant  and  very  slowly  varying 
errors  arising  from  products  of  oscillations  with  overlapping  frequency 
spectra.  Such  products  occur  as  unwanted  interactions  between  noise 
oscillations  and  stabilizing  coirpatations  or  errors  entirely  within  a 
particular  instrument.  Wideband  errors  such  as  gyro  drift  rate  grow 
proportional  to  the  square  root  of  time,  as  opposed  to  the  linear  growth  of 
constant  drift  rates. 

2.2  Global  Positioning  System  (GPS) 

The  Global  Positioning  System  (GPS)  is  a  system  of  navigational  satellites 
expected  to  reach  full  operational  3-D  capedsility  in  the  mid  1990s.  The  system 
configuration  will  consist  of  six  orbital  planes  inclined  55'  relative  to  the 
equator  at  eui  altitude  of  about  20,000  km  with  an  absolute  orbital  period  of  12 
hr.,  equivalent  to  a  revisit  time  to  the  seime  rotating  earth  position  of  24  hr. 

Each  orbital  plane  will  contain  three  active  satellites.  Three  additional 
satellites  will  serve  as  spares.  Usually,  at  least  four  satellites  will  be  in 
view  cinywhere  on  the  earth  at  any  time.  Presently,  seven  satellites  in  two 
orbital  planes  are  operational. 

It  is  expected  that  GPS  will  be  used  as  a  stand-alone  navigation  aid,  as  well 
as  with  INS,  heading  reference  systems,  altimeters,  on-board  precision  clocks, 
etc.  [Stein, 87]. 
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The  GPS  system  operates  by  broadcasting  four  types  of  information 
[Lundberg,87] ; 

o  time  signal 

o  pseudo  reinge  measurement 

o  integrated  doppler  measurement 

o  ephemeris  for  each  of  the  GPS  satellites 

Pseudo  reuige  is  obtained  by  converting  the  received  time-coded  navigation 
signals  into  offsets  relative  to  the  internal  clock  of  the  GPS  receiver.  The 
term  "pseudo  range"  is  used  because  of  the  synchronization  error  between  the 
user's  cind  the  GPS  clocks,  resulting  in  ambiguities  of  the  time  of  tremsmission 
of  the  received  signal. 

Two  L-band  signals  are  continuously  broadcast  by  each  satellite.  Two  coded 
versions  will  be  available  on  each  bemd,  but  the  higher  precision  P  code  will 
be  encrypted  whereas  the  lower  precision  C/A  code  will  be  generally  availadsle. 
The  pseudo  range  measurement  accuracies  will  be  on  the  order  of  1.5  m  (RMS)  for 
P  code  and  15  m  for  C/A  code  (Lundberg,87] ,  and  .1  nv/sec  velocity  accuracies 
for  P  code  [Teasley,87] .  The  integrated  doppler  measurement  accuracy  will  be  on 
the  order  of  .2  m  (RMS)  [Lundberg,87] . 

2.2.1  GPS  Error  Sources 

Four  types  of  degradation  are  the  major  error  sources  affecting  the 
determination  of  positional  accuracy  (Teasley,87] ; 

o  ionospheric  delay  errors 

o  tropospheric  delay  errors 

o  ephemeris  prediction  errors 

o  geometric  dilution  of  precision  (GDOP) 

Other  sources  of  error  include  (Kleusberg,87] : 

o  limited  measurement  resolution 

o  clock  errors 

o  receiver  noise 

o  signal  multipath 

Ephemeris  information  is  used  for  estimating  the  position  of  the  GPS  satellites 
in  view.  The  pseudo-range  and  integrated  doppler  measurements  are  used  to 
calculate  the  user's  position  relative  to  the  satellites  in  view. 
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The  accuracy  of  the  ephemeris  data  broadcast  by  GPS  is  claimed  to  be  about  12m 
(RMS)  if  ground  station  ephemeris  updates  are  regularly  received  three  times 
daily  by  GPS.  In  the  event  of  jamming  or  other  interference  with  these  updates, 
the  ephemeris  errors  could  grow  to  about  20  km  after  6  months  [Lundberg,87] . 
However,  by  using  cross-linking  of  range  aurid  integrated  doppler  signals  among 
the  GPS  satellites  themselves,  this  error  could  be  reduced  to  about  40m,  euid 
user  rcinge  errors  to  about  6m  [Menn,87]. 

Additional  loss  of  accuracy  can  occur  when  the  geometry  of  the  satellites 
visible  to  the  user  is  such  that  these  satellites  are  nearly  coplanar.  Another 
case  is  when  less  than  four  satellites  are  in  view  because  of  occlusion  by 
terrain  or  satellite  failure. 

Such  cases  of  unfavorable  geometry  are  often  quantified  by  the  "Geometric 
Dilution  of  Precision"  (GDOP).  GDOP  is  the  ratio  of  rms  position  and  clock 
error  to  the  1-sigma  satellite  rcinging  error  [Massatt,87] .  Another  ratio  which 
is  more  useful  to  most  users  is  the  "Positional  Dilution  of  Precision"  (PDOP). 
PDOP  is  the  ratio  of  rms  position  error  to  the  1-sigma  satellite  ranging  error. 
Other  ratios  for  particular  positional  component  errors  or  time  are  also 
defined.  A  more  intensive  discussion  of  these  issues  along  with  some  examples 
of  calculations  for  affected  latitudes,  longitudes  and  operating  times  given 
minimum  accept2±)le  values  of  PDOP  is  found  in  [Massatt,87 ) . 

The  first  case  occurs  for  a  few  regions  along  the  40th  parallel.  Some  of  these 
situations  could  be  improved  using  equatorial  orbits  for  the  three  spares 
[Stein, 87] . 

Models  of  tropospheric  delay  and  dual  frequency  measurements  may  be  used  to 
reduce  the  ionospheric  delay  errors  (Kleusberg,87] . 

2.2.2  Differential  GPS 

There  are  a  number  of  methods  for  potentially  increasing  the  accuracy  of 
positional  determination  using  either  the  doppler  portion  of  the  signals  or 
positional  data  relative  to  a  known  position.  Such  methods  are  referred  to  as 
"Differential"  GPS. 

One  such  scheme,  described  in  [Beser,87],  proposes  the  use  of  a  stationary  GPS 
receiver  at  a  previously  surveyed  monitor  location.  Using  the  GPS-computed 
position  and  the  known  position  of  this  monitor,  correctional  estimates  can  be 
computed  and  broadcast  to  users  by  GPS. 

Such  a  relatively  simple  scheme  would  work  reasonably  well  for  three  reasons 
[Kleusberg,87] .  First,  the  satellite  clock  errors  for  a  remote  receiver  and 
monitor  would  be  the  same  and  are  cancelled.  Also,  satellite  ephemeris  errors 
affect  measurements  of  simultaneous  observers  almost  equally.  Therefore,  the 
computed  position  of  a  remote  receiver  could  be  computed  relative  to  the 
monitor  position.  Finally,  the  remaining  errors  due  to  receiver  noise,  limited 
measurement  resolution,  and  signal  multipaths,  are  more  in  the  high  frequency 
range.  These  high  frequency  errors  can  be  reduced  by  either  smoothing  with  INS 
data  or  using  the  GPS  carrier  signal  doppler. 
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Two  alternatives  for  the  use  of  doppler  data  are  presently  considered.  One  such 
method  would  filter  the  pseudo  range  time  series  with  the  doppler  measurements 
for  each  visible  satellite  in  order  to  derive  the  receiver's  position  euid  clock 
offsets.  Another  alternative  would  use  the  doppler  measurements  to  filter  the 
receiver  positions  which  have  been  separately  derived  from  pseudo  range 
measurements . 

Both  of  these  methods  would  work  better  if  smoothing  of  the  entire  appropriate 
data  set  were  employed  rather  than  recursive  filtering. 

2.3  Combined  Use  of  INS/GPS  for  Control  Information 

A  short  summary  comparison  of  the  relative  advantages  of  INS  vs.  GPS  is  given 
in  Table  2.1  [Teasley,87 1 .  It  is  clear  that  there  would  be  advantages  in  using 
the  two  data  sources  together  to  reduce  the  effect  of  INS  drift  errors,  GPS 
jamming,  and  GPS  tracking  of  high  dynamic  maneuvers.  Other  uses  include 
initialization  of  the  INS  alignments  using  GPS,  using  INS  in  the  presence  of 
partial  or  degraded  GPS  signals,  and  using  the  INS  to  aim  the  GPS  antenna  and 
control  the  receiver  bandwidth  [Bose, 87]. 

The  types  of  GPS/INS  combinations  are  discussed  in  section  2.3.1,  a  scheme  for 
loosely  coupled,  flight  path  estimation  using  combined  GPS/INS  data  appears  in 
section  2.3.2,  and  an  example  of  the  use  of  this  procedure  on  real  data  is 
given  in  section  2.3.3.  A  brief  discussion  of  some  of  the  issues  for  real-time 
compensations  to  raw,  complex  SAR  data  is  given  in  section  2.3.4. 

2.3.1  Types  of  Combined  GPS/INS 

INS  and  GPS  data  can  be  combined  jointly  using  a  number  of  processing  options. 
These  options  can  be  classified  as  loosely  coupled  or  strongly  coupled 
[Teasley,87 ] . 

Loosely  coupled  processing  essentially  involves  processing  the  GPS  data 
separately,  and  then  using  it  to  compensate  for  the  bias  arising  from  the  INS 
drift  errors.  Therefore,  a  Kalman  filter  formulation  for  processing  the  INS 
data  measurements  models  augmented  states  which  include  bias  terms. 

However,  there  are  problems  arising  from  the  different  clock  errors  for  the  two 
time  sources. 

Tightly  coupled  processing  involves  processing  the  raw  GPS  and  INS  data  in  a 
single  filtering  operation.  The  GPS  pseudo  range  and  doppler  measurements  and 
the  INS  linear  and  angular  acceleration  measurements  are  filtered  using  a 
combined  dynamic  model.  In  addition,  the  short-term  INS  data  can  be  used  to 
track  the  GPS  carrier  signals. 

A  simulation  of  such  combined  processing  is  reported  in  [Bose, 87).  The  initial 
error  sources  and  some  representative  error  budgets  are  given  in  Tables  2.2, 
2.3,  2.4,  and  2.5.  Results  are  given  for  the  non- jamming  case  in  Table  2.6  and 
for  the  case  of  GPS  jamming  in  Table  2.7. 
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Interestingly,  one  conclusion  of  this  study  [Bose, 87]  was  that  strapdown 
systems  were  found  to  be  superior  in  performance  for  land  vehicle  navigation. 
This  was  mostly  because  the  z-gyro  height  bias  was  more  observable  over  hilly 
terrain  using  strapdown  rather  than  gimballed  systems. 

2.3.2  Flight  Path  Estimation 

This  section  contains  an  estimation  procedure  for  combining  GPS  emd  INS  in  a 
loosely  coupled  form  to  estimate  a  flight  path.  This  procedure  is  based  on 
Kalmzm  smoothing,  as  distinct  from  Kalman  filtering.  Essentially,  such  a 
smoothing  procedure  allows  all  the  GPS  and  INS  data  for  a  given  imaging  period 
to  be  used  as  a  whole  rather  them  only  successive  portions  of  it  to  be  used 
recursively.  That  such  processing  should  lead  to  better  estimates  is  not  only 
intuitively  appealing,  but  also  theoretically  valid  [Bierman,73] ,  [Gelb,741, 
[Maybeck,79] . 

Although  such  a  procedure  cannot  be  as  efficiently  implemented  "on  the  fly"  as 
a  Kalman  filter,  it  can  be  computed  as  a  combination  of  time-forward  and  time- 
backward  Kalman  filters.  Of  course,  all  of  the  processed  time-forward  data  must 
be  kept  to  implement  the  time-backward  filter. 

The  following  is  an  implementation  of  a  smoothing  procedure  based  on  the 
computational  approach  of  [Bierman,731  for  using  GPS  and  INS  data  to  estimate 
flight  path  parameters. 


Algorithm  to  estimate  aircraft  position  from  INS/GPS  data 
1 .  Coordinate  System 

Let  R  =  (Rj, ,  Rj. ,  )  be  the  component  of  the  aircraft  position  vector  as 

measured  in  a  rotating  geographic  coordinate  frame  with  orthogonal  unit 
vectors  I^.  ,  ,  K^ .  The  origin  is  fixed  to  the  earth  at  latitude, 

longitude  (X,<]>);  I^  points  towards  local  North;  J^,  points  toward  local 
East;  Kg  points  down. 

The  apparent  aircraft  velocity  as  measured  in  this  same  frame  is  V  =  (V^, , 
Vg,  V^);  i.e.,  Vj,(t)  =  d  Pf,(t),  etc 

at 


2.  Single  channel  dynamic  model. 

Consider  the  case  of  1  accelerometer  which  measures  acceleration  in  the  Ig 
(North)  direction.  Let  the  state  vector  be 


X  = 

LvJ 

For  short  time  durations  we  use  the  following  model  [Farrell,  76]; 
(1)  X{t)  =  FX(t)  +  u(t) 
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where  F  is  a  (2x2)  matrix  of  coefficients  and  u  is  a  (2x1)  vector 
representing  the  forcing  fianction.  Since: 

^  (t)  -  V„  (t) 

V(,  (t)  »  Ajj  (t)  (the  measured  acceleration), 

we  have 


Equation  (1)  has  the  solution 

x(t)  -  ♦(t,tg)x(tg)  +  *(t,T)  u(T)  dr 

where  the  transition  matrix  4  is  computed  by 

♦  (t,T)  =■  e''<'=-'f' 

=-  I  +  F*(t-T)  +  F^  (t-T)2  +  ... 

-  I  +  F  .  (t-T)+0 


23 


Since  the  actual  motion  satisfies:  X  -  Fx  +  u 

and  the  measured  motion  satisfies:  X  -  Fx  +  u 

the  errors  must  also  satisfy:  X  -  Fx  +  w 

At  times  {t^,  k  -  1,..N}  we  get  measurements  of  position  from  GPS 

where  (t^)  is  the  true  position,  eund  V,,  is  the  error  in  the  GPS  position 
measurement.  At  each  time  t^ ,  form  the  difference  between  the  INS 
measurements  auid  the  GPS  measurements: 

Zk  -  )  -  g, 

-  (Rr,(tk)  -  Rn  (tj)  -(R^(t,)-V,) 

“  ^  ^  ^k  ) 


where  we  have  now  put  this  in  the  Kalman  filter  format.  The  "measurement 
matrix"  is  H  -  (-1  o].  The  "measurement  error"  is  the  error  in  the 

GPS  positioning.  It  is  assumed  that  E(V^ )  -  0  and  CoVar  (V^ )  -  . 

The  system  model  is 

x(t)  -  Fx(t)  +  w(t) 
or  p^(t)1  To  11  p^(t)l  To 

.V„(t)J  Lo  oj  LV„(t)J  L\(t). 

Aj,(t)  is  noise  in  the  INS  acceleration  measurements,  with 
E(A„(t))-0,  Var  (A„(t))  -  . 

The  covariance  of  w(t)  is  therefore 

o  o 

Q( t)-covar(w( t) )  - 

.0 
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4.  Single  channel  estimation  procedure. 

Since  all  of  the  measurements  needed  by  the  estimation  process  are 
available  before  beginning  the  procedure,  we  can  use  a  Kalmaun  smoother. 
This  is  guarainteed  to  give  smaller  errors  them  the  "real  time"  Kalnan 
filter  [Gelb,  74],  [Maybeck,  82).  The  overall  plain  to  estimate  the 
aircraft  position  is  as  follows: 

a)  Use  the  GPS  position  estimates,  assumed  GPS  measurement  variances,  eund 
assumed  INS  accelerometer  variance  to  estimate  the  values  of  the  INS 
error  vector 

x(t).  Use  the  Kalmain  smoother  in  step  5  below. 

b)  Form  the  estimate  of  position  and  velocity  using 

x(t)  -  x(t)  +  x(t) 

where  x(t)  is  the  measured  INS  data.  The  covariance  of  these  estimates  is 
given  by  the  U^  matrices  in  the  Kalman  smoother. 

5.  Kalmeui  Smoother  (Gelb  74],  [Maybeck,  82] 

System  Model 


x(t)  -  F  x(t)  +  w(t) 


where  x 


'0  1' 
.o  o. 


covar  (w( t) )- 


o  o 
Lo  T.V 


Measurement  Model 

Zk  ”  Hk^k  + 

where  -  [-1  o] 
w^  =-  x(tk  ) 
var(v,  )  - 

Initial  conditions 


E(x(t,  )  )  =  X,  , 

covar  (x(  t,  ) )  =  =  P,' 

E(w(  t^  )v^ )=o 

for  k  »  1  to  N,  begin  (iteration  on  k  for  forward  filtering  step) 


25 


state  estimate  extrapolation 


In  the  following,  the  notation  (-)  refers  to  a  quantity  before  am  update,  and 
(+)  after  an  update. 

States  are  first  extrapolated  from  the  (k-1)^  stage  to  the  stage,  and  then 
estimated  at  the  k^  stage. 


♦k-i  J'k-i  (  +  )  ^ere 


^  ^-k  +  l  ^ 

.O  1  . 


Pk(-)  -  ♦k-i  Pk-1  <+)  ♦k-i  +  Qk 

Kk  -  P,(-)  [H,P,(-)H,t+R^]-i 


State  estimate  update 


Kk  -  Pk(-)Hx"  [H,P,(-) 

x^(t)  -  X^(-)  +  (z^-H^x^(-)] 

Pk(  +  )  -  [I-K^HJP^(-) 

For  backward  filtering  step, 
set  y^-x^" 

Un  - 

for  k  »  N-1  down  to  1,  do  begin 
A,  -  PkM  +  )^^P, 

Yk  -  Xk(  +  )  +  \  tYk^i  - 
Uk  -  Pk(+)  +  \  tUk.i  -  Pk*i 


end 

{y^ ,  k  *  1,...N}  are  the  resulting  smoothed  estimates 
{U^ ,  k  •  1,...N}  are  the  smoothed  covariances 

6 .  3-dimensional  extension 

The  simple  model  given  here  may  be  generalized  to  each  of  the  three 
dimensions.  The  resulting  dynamic  model  is  uncoupled.  If  one  also  assumes 
that  the  measurement  errors  on  each  of  the  three  directions  I^ ,  are 

uncorrelated,  one  can  solve  for  each  axis  individually. 
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2.3.3  Example  Data  Set 


An  example  of  the  use  of  the  smoothing  method  of  the  previous  section  may  be 
demonstrated  on  Intera's  STAR2  imagery  in  a  subsequent  Option.  Both  INS  2uid  GPS 
data  are  available  with  the  processed  SAR  imagery.  The  INS  and  GPS  data  were 
not  filtered  together  by  Intera  using  a  tightly  coupled  algorithm,  but  were 
processed  separately.  We  will  combine  the  two  data  sources  by  smoothing  in  a 
loosely  coupled  fashion. 

The  INS  data  was  used  to  make  real-time  phase  and  rauige  corrections  to  the  raw 
SAR  data  in  the  SAR  processor  (see  section  2.3.4). 

2.3.4  Real-Time  Compensations  for  SAR  Imaging 

Real-time  compensations  for  aircraft  SAR  imagery  differ  from  those  used  in 
satellite  SAR  imagery  (see  sec.  5.)  because  of  the  differences  in  platform 
stability,  the  importcince  of  earth  rotation  effects,  and  other  factors.  In 
Option  1,  these  effects  may  be  studied  in  greater  detail. 
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COMPARISON  OF  INERTIAL  NAVIGITION  SYSTEM 
AND  GLOBAL  POSITIONING  SYSTEM 


Parameter  GPS 


Position  X 

Velocity  X 

Time  X 

Acceleration 
Attitude 

High  Dynamic  Performance 
Time  Invariant  Accuracy  X 


Immunity  to  Countermeasures 


Table  2.1  (from  (Tensley,  87)) 


GPS  EPROR  SOURCE  BUDGET 


GFS  Error  Source 

Syr.bol 

SufTtUe  r 
of 

Sources 

SPACE  VErilC'.E 

Clock 

Phase  Ir.ltlal  condUron 

'  'SO 

U 

prequenoy  bias 

;fst> 

It 

Frequency  -■hire  noise 

'  *  sw* 

** 

Ephenerls 

12 

position  ir.ltlal  ccnbltlon 

i  y.gQ  s  1  ^  so  • '  “SO 

?csl:icn  rate  bias 

-'■srb  '-^  srb  ‘•srb 

;  2 

ATMCSrrLrlC 

Trepesphere 

Phase  bias 

-  ’tb 

Correlated  noise 

'  *  t  c 

loncspr.e  re 

Phase  correlated  noise 

'  'l 

hultlqach 

Foslrlon  white  noise  sar.qles 

etnu 

USER 

Clock 

Phase  initial  condition 

^  '-JO 

1 

Frequency  initial  condition 

fjuo 

1 

Frequency  rate  bias 

-  -  rub 

‘ 

Frecjcncy  white  noise 

-  •  uw 

Frequency  rate  white  noise 

*  •  r  'jw 

1 

receiver 

“ j  1 1 1  pa th/pseuderar ce  white  noise  sarples 

■  -■•pw 

Celta  pseudcTar.ie  ■-‘‘.ite  nclse 

Table  2.2  (from  (Bose,  87]) 
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GYROSCOPE  ERROR  BUDGETS 


■■ 

■H 

RMS  \'alues 

Lo»' 

Oiialily 

Medium 

Quality 

High 

Quality 

Rian 

'’,Hr 

1.0 

0.1 

0.01 

Scale 

Faciof 

rr 

0.1 

0  02 

0.005 

Scale  Factor 
Asymmetry 

/f 

001 

0.002 

0.0005 

Mass 

Unbalance 

°'nr;g 

0.5 

0.1 

0.02 

Quadrature 

°'Hr'g 

0.3 

0  05 

0  01 

j  Aniso- 
elaslicity 

0.: 

0  04 

0  01 

1 

1  Nonorthc- 
!  gonality 

sTc 

100 

20 

10 

1  Correlated 
i  ’^cine 

HfS  ffc 
(CT  .1 

0.03  <2  300 

0  01  '«  200 

0  002  'S  200 

1  White  N'oiie 

1 

‘^'Hr  -/hT 

0.: 

0.1 

0  05 

Trend 

“'Hr 'Day 

0  01 

0  002 

0.0005 

Table  2.3  (from  (Bose,  87]) 


ACCELEROMETER  ERROR  BUDGETS 


RMS  N'alufS 

Error 

S'^urcj 

Units 

Low 

Quality 

Medium 

Quahry 

H^ch 

0',:s:ily 

Bias 

ue 

.'00 

100 

50 

Scale 

Factor 

0  2 

0  V-1 

0  01 

Scale  Factor 
Asymmetry 

0  05 

0  01 

0  CO  5 

Direct 

Quadratic 

Nonlinearity 

C- 

140 

‘0 

20 

Cros! 

Quadralic 

Noniincirity 

;.■£  r 

50 

20 

< 

Nonorthc- 

gcnality 

sec 

■jO 

:o 

:0 

Ccnelaied 

Norse 

^  C  T 

10'^  20  se  r 

5  20  sec 

1  20  '?c 

i 

hiie 

N'oise 

-E'  vHi 

-•0  '  ! 

■0 

?  1 

1 

Trtnd 

ug  Day 

5  ! 

“  1 

•  ! 

1 

1  I 

Table  2.4  (from  (Bose,  87]) 
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INERTIAL  SYSTEM  INITIAL  CONDITION  AND 
GRAVITY  DISTURBANCE  ERROR  BUDGETS 


Error 

RMS 

Soute* 

Symbola 

Uniti 

Values 

Level  Pos) tion 

‘fyo 

Sec 

3  0 

Altitude 

thQ 

Feet 

0 

Velocity 

iV.O-  ''vO- 
■ 

Ft/sec 

0.1 

VerticAlity 

'.O'  'yO 

Decree  a 

1.0 

Azimuth 

*l0 

Deoreet 

2.0 

Gravity 

Deflection 

>1,.  «Sy 

9  NH  (C. C. ] 

35  9  20 

Gravity  Anomaly 

vq  f  NM  <C. D. 1 

3  5  9  20 

Deflection  Initial 
Condition 

L-9 

' 

Ancrraly  initial 
Condition 

-9 

35 

Table  2.5  (from  [Bose,  87)] 


PERFORMANCE  SUMMARY  OF  LOW,  MEDIUM,  AND  HIGH  QUALITY 
HYBRID  AIDED  INERTIAL  NAVIGATION  SYSTEM  CANDIDATES 


'-'vbrid  Svsiem 
CandiCatei 

Average  RSS 
East  Velocity 
Ef'or 

(Mett!i,'Sec\ 

Ave'age  RSS 
North  Velocity 
Error 

(W.tttts'Setl 

Average  RSS 
Vertical  Velocity 
Error 

(Wtie'y'SecI 

Average  RSS 

Ajipnuth 

Error  [Mils) 

Peak  RSS 
Horitonial  Position 
CEP  IMetersI 

Peak  RSS 
Attitude 
Error  (Meters) 

GIMBALLEO  INERTIA 

L 

t  Low 

0,120 

0.122 

0.037 

21.8 

360.0 

38.7 

Odo.Tieter  Only 

[  Medium 

0.049 

0.C55 

0.030 

2  87 

62. 0 

38.7 

High 

0.046 

0.049 

0027 

0,74 

72.0 

38.0 

/  Low 

0.040 

0  043 

0  037 

100 

52.5 

38.7 

Wild  PLRS 

Vedium 

0,03 

0.C37 

003 

1,T5 

25.8 

38.7 

[  High 

0,03 

0,037 

0  027 

1  20 

26.5 

38.0 

/  Low 

0  012 

0  015 

0015 

4  9 

6  0 

13.2 

W.th  GPS 

Medium 

0.01 

0,012 

0015 

0.79 

6,0 

13.2 

High 

0,0091 

0-01-3 

0  0146 

0.46 

6  0 

13  2 

f  Low 

0  012 

0  015 

0.015 

4  9 

5  5 

13.2 

With  GPS  S  PLRS 

Medium 

0,01 

0.012 

0015 

0,79 

5  5 

13,2 

i  High 

0.0091 

0  0113 

00146 

0  46 

5.5 

13.2 

strapdown  inertial 

c  Low 

0.076 

0.091 

0  040 

10.2 

131  0 

32.3 

Odometer  Only 

Medium 

0.046 

0.046 

C030 

1,53 

66  0 

31.7 

High 

0.043 

0,043 

OC30 

0  69 

62  0 

31.7 

Low 

0.040 

0  046 

0  037 

7  64 

41  6 

32.3 

With  PLRS 

Med»um 

0.03 

0  037 

0.C3 

1  30 

20.3 

31.7 

High 

0.03 

0  037 

CC3 

0.74 

17.7 

31.7 

{ 

Low 

0  015 

0  017 

0015 

5.1 

6.1 

13.7 

With  GPS  ( 

Medium 

0.01 

0  012 

0.015 

0  83 

5.0 

13.7 

1 

High 

0  00S5 

0  0116 

0.015 

C.44 

6  0 

13.7 

Low 

0.015 

0,017 

0015 

5.1 

5  5 

•3.7 

With  GPS  i  PLRS 

Medium 

0.01 

0  012 

00-5 

0  33 

5  5 

■3,7 

Htgh 

0  0055 

0  01-6 

0015 

C  44 

5  5 

13  7 

Table  2,6  (from  [Bose,  87]) 
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HYBRID  SYSTEM  PERFORMAMCE  DEGRADATION 
DUE  TO  ONE  HOUR  TO  .L  JAMMING 


Peak  RSS 

Peak  RSS 

Average  RSS 

Hybrid  Syitem 

Horizontal  Position 

Altitude 

Azimuth 

Candidates 

CEP  (Meters) 

Error  (Meters) 

Error  (Mils) 

Unjammed 

Jammed 

Unjammed 

Jammed 

Unjammed 

Jammed 

GIMBALLED  INERTIAL 

Low 

52.5 

68.0 

38.7 

■I 

9.95 

12.04 

With  PLRS 

f'/.edium 

25.8 

47.0 

38.7 

Km 

1.76 

2.41 

High 

26.5 

50.0 

38.0 

1.20 

1.90 

f  Low 

6.0 

140.0 

13.2 

31.2 

4.86 

10.18 

With  GPS 

Medium 

6.0 

36.5 

13.2 

31.1 

0.79 

1.81 

High 

6.0 

33.0 

13  2 

31.1 

0.46 

0.97 

f  Low 

6.0 

16.0 

13.2 

31.1 

4.86 

8.33 

With  GPS  8< 

Medium 

6.0 

8.8 

13.2 

31.1 

0.79 

1.44 

Unjammed  PLRS 

High 

6.0 

8.8 

13.2 

31.1 

0.46 

0.E3 

STRAPDOWN  iNERTlAL 

Law 

41.6 

53.9 

32.3 

32.6 

K 

8.56 

With  PLRS 

Medium 

20.3 

37.0 

31  7 

31  8 

1.39 

High 

17.7 

33.6 

31.7 

31.3 

0.74 

0.74 

Low 

6  0 

75.0 

13.7 

20  8 

5.09 

7.64 

With  GPS 

Medium 

6  0 

36  5 

13.7 

20.7 

0  83 

1.25 

High 

6  0 

36.5 

13.7 

20.7 

0  44 

0.63 

With  G=‘S  &  j 

Low 

5  5 

14.0 

20  8 

5  09 

6.48 

Unjammea  PLRS 

’/ed'j.Ti 

High 

5  5 

5  5 

8.1 

8.1 

_ 

20.7 

20.7 

0  83 

0.44 

1.11 

0,60 

Table  2.7  (from  [Bose,  87]) 


3.0  Generation  of  Control  Information  Using  Terrain  Data 


This  section  discusses  the  use  of  terrain  data  for  extracting  control 
information  to  calibrate  imagery.  The  simulated  scenario  will  assume  that 
approximate  position  estimates  of  the  sensor  platform  will  be  available  from  an 
INS  or  combined  INS/GPS  source,  and  that  the  flight  path  direction  will  be 
accurate  to  within  a  degree  [Mercer, 87]. 

Two  distinct  methods  for  extracting  control  from  terrain  data  will  be  described 
here.  Both  methods  will  employ  synthetic  images  generated  from  assi:iined 
flightlines  using  terrain  models. 

Because  the  available  sensed  images  are  ground  range  images,  the  synthetic 
images  must  also  be  ground  range  images  for  emy  comparisons  to  occur.  Because 
sleint  range  images  conpress  the  near-ranges  compared  to  the  far  ranges,  ground 
range  images  are  produced  vdiich  are  more  visually  pleasing.  However,  there  are 
errors  associated  in  producing  the  conversion  to  ground-range,  and  these  errors 
are  dependent  on  the  terrain.  Therefore,  accurate  ground  range  presentation 
really  requires  registration  of  the  slant  range  image  with  the  terrain. 

The  approach  for  generating  control  information  for  a  SAR  image,  i.e.  assigning 
a  geographic  location  to  every  pixel,  will  be  to  determine  the  sensor  platform 
position  as  a  function  of  time.  Ideally,  the  flight  path  for  aui  aircraft  should 
be  straight  and  level.  Of  course,  this  does  not  occur  in  practice  (see  section 
2.3.3  and  2.3.4).  Therefore,  straight  line  flight  path  approximations  for 
smaller  blocks  of  azimuthal  lines  will  be  attempted.  These  straight  line 
approximations  can  be  used  as  supplementary  update  measurements  to  remove  the 
bias  effects  of  drift  in  the  INS  data  for  estimating  a  flight  path  as  a 
function  of  time.  In  this  way,  these  measurements  could  be  used  in  the  same  way 
as  the  GPS  measurements  are  used  for  updates  in  the  Kalman  smoothing  algorithm 
described  in  section  2.3.2.  Of  course,  GPS  measurements  would  be  much  more 
accurate. 

Once  a  flight  path  or  smaller  flight  path  segments  have  been  determined,  i.e. 
the  "resection  in  space"  has  been  estimated,  the  range  and  doppler  equations 
for  each  pixel  can  be  used  to  intersect  with  the  terrain  model  to  obtain 
geographic  coordinates  [Leberl,83],  [Kwok,  87],  [Curlander,  87).  In  the  case  of 
airplane  SAR,  usually  no  squint  is  used,  and  the  doppler  cone  degenerates  into 
a  sccinning  plane. 

An  alternate  approach  would  be  to  try  to  register  the  image  to  model 
coordinates  without  determining  the  exterior  orientation  of  the  sensor,  i.e. 
the  flight  path.  As  is  discussed  in  [Leberl,83],  such  approaches  generally  fall 
into  three  categories: 

o  "riibber  sheeting"  interpolative  methods  using  ground  control  points 

o  using  the  projection  equations  and  ground  control  points  to  estimate 
parameters 

o  hybrid  combinations  of  the  above 
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It  would  seem  that  using  the  model  to  obtain  the  sensor  flight  path  is  the 
better  approach  because  it  directly  uses  the  model  information  to  generate 
control  features  such  as  shadows  or  shading. 

The  generation  of  synthetic  SAR  images  appears  in  section  3.1.  The  use  of 
shadows  created  by  terrain  occlusions  for  estimating  sensor  platform  position 
is  discussed  in  section  3.2.  The  use  of  terrain-induced  shading  for  estimating 
sensor  platform  position  is  discussed  in  section  3.3.  Finally,  a  brief  summary 
of  the  problem  of  resection  in  space  is  treated  in  the  Appendices,  sections  9.1 
and  9.2. 

3.1  Generation  of  Synthetic  SAR  Imagery 

Synthetic  image  generation  for  both  SAR  and  optical  imagery  is  discussed  in 
this  section. 

3.1.1  Synthetic  SAR  Imagery 

The  present  software  for  forming  projected  views  assumes  that  the  flight  line, 
and  therefore  the  azimuthal  direction,  is  parallel  to  one  of  the  model 
coordinate  euces.  If  this  is  not  the  case,  then  the  appropriate  rotation  of  the 
model  must  be  performed.  Operationally,  this  is  not  much  of  a  constraint  since 
the  flight  path  is  assumed  to  be  known  within  a  degree.  No  squint  angle  is 
ass\amed,  which  is  the  usual  case  for  aircraft  SAR  imagery.  Therefore,  the 
doppler  equations  reduce  to  planar  sccinning. 

Essentially,  the  formation  of  projected  views  requires  the  calculation  of  a 
local  surface  orientation  of  the  terrain  model  for  each  resolution  element,  and 
the  suppression  of  occluded  regions.  Another  component  of  image  intensity  which 
affects  gray  values  in  a  real  sensed  image  is  local  thematic  value,  ie.  the 
scattering  cross  section  coefficient. 

Because  such  thematic  values  are  not  available  in  most  terrain  model  data, 
their  effect  is  ignored  in  these  simulated  views.  The  use  of  thematic  values  is 
irrelevant  for  the  calculation  of  shadows  in  synthetic  images. 

However,  when  non-shadowed  regions  in  a  synthetic  image  are  to  be  compared  to 
the  corresponding  regions  in  a  real  sensed  image,  there  will  generally  be 
intensity  distortions  which  are  not  related  to  terrain.  This  issue  is  discussed 
in  section  3.3. 

The  geometry  of  the  imaging  scenario  and  its  inputs  is  shown  in  Fig.  3.1. 

There  are  a  number  of  considerations  for  producing  such  synthetic,  aircraft  SAR 
images: 

o  intensity  values 

o  occlusions/radar  shadows 

o  range  overlay 

o  discretization  artifacts 
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The  intensity  values  are  computed  usii  ^  the  cosine  of  the  local  incidence  angle 
of  the  radar  beam  with  respect  to  the  local  surface  facet  normal.  These 
intensities  are  computed  for  every  terrain  cell  surface  facet. 

A  separate  file  for  occlusions  or  overlay  is  computed  for  every  terrain  cell 
surface  facet.  Occlusions  are  computed  simply  by  checking  for  changes  in  the 
monotonicity  of  the  angle  (see  Fig.  3.1)  as  one  proceeds  from  the  far  range, 
rj ,  to  the  near  range  r  .  Overlay  is  determined  by  whether  the  monotonicity  of 
reunge  offsets  is  violated  as  one  proceeds  from  far  to  near  reuige. 

Discretization  artifacts  can  potentially  occur  in  the  construction  of  a 
synthetic  image  because  the  changes  in  the  range  for  terrain  cell  increments 
may  involve  jump  discontinuities  when  measured  as  range  resolution  increments. 
If  uncorrected,  this  leads  to  a  synthetic  image  for  vdiich  there  are  no  returns 
for  certain  range  resolution  pixels. 

A  simple  averaging  scheme  using  oversampling  is  employed  to  offset  this 
anomaly.  Basically,  this  procedure  oversamples  the  range  buffer  by  some  factor. 
For  practical  purposes,  eui  oversanpling  factor  of  4  was  sufficient.  Then  the 
intensities  for  the  final  range  cells  of  correct  resolution  are  computed  by 
averaging  groups,  in  this  case  groups  of  4,  of  contiguous  cells  from  the 
oversampled  buffer. 

3.1.2  Synthetic  Optical  Imagery 

A  freuned  optical  image  is  formed  by  pixels  representing  azimuthal  and  elevation 
angle  "rays",  instead  of  pixels  representing  azimuthal  and  range  distances  as 
in  SAR.  Ttie  same  cosine  relationship  is  used  for  calculating  intensities  as 
above.  No  overlay  occurs  as  in  SAR  imaging,  but  shadows  and  occlusions  are 
calculated  using  monotonicity  of  projection  angle  as  before. 

3.2  Use  of  Terrain-Induced  Shadows  for  Control 

There  is  no  problem  determining  the  occurrence  of  a  shadow  or  locating  its 
boundary  in  a  synthetic  image.  However,  the  situation  is  obviously  not  quite  as 
straight-forward  in  real  imagery. 

Shadows  in  a  real  SAR  image  are  subject  to  two  kinds  of  instabilities.  A 
geometrically-induced  instability  involves  the  problem  of  edge  migration  due  to 
small  perturbations  in  the  platform  flightpath  and  unmodeled  perturbations  in 
the  real  terrain.  A  radiometrically-induced  instability  involves  the 
uncertainty  in  the  choice  of  a  precise  shadow  threshold  because  of  the  effects 
of  noise  cind  speckle. 

These  geometrically-induced  instabilities  can  arise  whenever  there  is  a  portion 
of  real  terrain  present  that  is  not  convex.  Departures  from  convexity  can 
potentially  introduce  relatively  large  changes  in  shadow  contours  for 
relatively  small  chcmges  in  platform  position. 


The  extent  of  shadow  boundary  instabilities  due  to  changes  in  the  shadow 
threshold  are  potentially  greater  in  the  side  of  the  shadow  furthest  away  from 
the  sensor.  The  shadow  side  closest  to  the  sensor  should  have  a  stronger 


radiometric  edge  because  the  terrain  slopes  near  this  geometric  discontinuity 
often  approach  a  near-specular  orientation. 

Because  of  the  potential  for  such  large  changes,  precise  measurements  on 
shadows  are  more  appropriate  for  finer  estimation  of  platform  position  rather 
them  for  acquisition. 

3.2.1  Extraction  of  Shadows  From  Imagery 

The  approach  for  extracting  the  shadowed  regions  in  the  real  imagery  involves 
thresholding  and  vector ization.  These  issues  are  described  in  sections  3. 2. 1.1 
and  3. 2. 1.2  respectively. 

3. 2. 1.1  Thresholding 

Segmenting  shadows  involves  discriminating  between  highly  specular  no-return 
areas  such  as  large  bodies  of  water  and  true  radar  shadows  due  to  terrain 
occlusions.  It  also  involves  carefully  choosing  a  shadow  threshold  so  that  a 
clean  delineation  of  the  shadow  boundary  from  the  surrounding  region  is  formed. 

Solving  the  former  problem  of  determining  a  specular  no-return  region  is 
greatly  aided  by  the  use  of  prior  knowledge  of  the  nearby  existence  of  lakes, 
reservoirs,  etc.  Such  knowledge  would  be  represented  in  DFAD  or  other  feature 
data.  This  feature  would  then  be  actively  searched  for  in  the  saime  way  that  a 
particular  terrain  shadow  is  searched  for,  given  some  hypothesized  sensor 
platform  position. 

Without  such  prior  knowledge,  an  evaluation  of  the  degree  of  "goodness  of  fit" 
of  each  segmented  potential  shadow  region  would  be  the  only  way  of  deciding 
that  the  segmented  region  is  a  no-return  region,  or  that  there  are  terrain 
model  errors.  Such  decision  issues  can  be  approached  using  a  rule-based 
approach,  and  are  discussed  in  section  7. 

An  excimination  of  the  intensities  of  a  real  STARl  Intera  SAR  image  (see  Fig. 
3.4)  shows  a  clear  bimodality  of  the  histogram  for  shadow  and  non-shadow 
regions.  This  is  shown  in  Fig.  3.5.  The  variations  in  shadow  outline  are  shown 
in  Fig. 3. 6  -  3.10.  As  can  be  seen  in  these  images,  the  variations  in  shadow 
shape  are  reasonably  well-conditioned  with  respect  to  small  threshold  changes 
around  the  local  minimum  in  the  histogram  separating  shadow  from  non-shadow 
intensities. 

For  SAR  images  with  a  great  deal  of  noise  and  speckle  content,  it  may  be 
necessary  to  prefilter  the  images  using  some  nonlinear  spatial  filter  prior  to 
calculating  the  histogram.  Two  examples  of  such  nonlinear  filters  which  reduce 
these  effects  but  still  preserve  edge  structures  relatively  well  are  described 
in  [Lee, 86].  One  of  these,  the  "sigma  filter",  has  been  used  by  VEXCEL  on 
SEASAT  SAR  images  of  arctic  ice  floes  to  reduce  the  effects  of  speckle  prior  to 
thresholding  using  the  image  histogram. 

The  sigma  filter  is  based  on  a  Gaussian  distribution  model  of  the  pixel 
intensities  in  a  (5x5)  or  (7x7)  pixel  window.  The  algorithm  uses  this 
statistical  model  for  selective  averaging  of  pixel  intensities.  The  assumption 
is  that  real  contours  arise  as  pixels  representing  the  higher  percentiles  of 
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the  local  distribution.  Unfortunately,  spot  noise  also  satisfies  this 
criterion.  Only  pixels  belonging  to  the  two-sigma  range  of  the  normal 
distribution,  i.e.  99.5%,  are  averaged.  The  pixels  falling  outside  this  range 
remain  as  before.  A  special  thresholding  procedure,  somevdiat  ad-hoc,  is  then 
applied  to  remove  spot  noise. 

The  other  of  these  two  algorithms  is  called  the  "local  statistics"  method. 
However,  this  neime  is  somewhat  misleading,  since  the  real  basis  of  this 
algorithm  is  estimation  theory  using  local  statistics  in  (5x5)  auid  (7x7) 
neighborhoods.  The  method  is  based  on  a  multiplicative  noise  model  relating  the 
observed  pixel  intensity  as  the  product  of  the  true  "signal"  at  that  pixel 
location  and  eui  independent  noise  term.  A  linear  series  exp2msion  of  this  model 
giving  the  observed  intensity  as  a  function  of  the  meam  values  of  both  the  true 
signal  and  the  noise  term,  assumed  to  be  unity,  allows  an  estimation  of  the 
meein  and  variance  of  the  true  signal  as  a  function  of  the  mean  and  variance  of 
the  observed  intensity.  These  estimates  of  the  mean  and  variance  of  the  true 
signal  are  then  used  in  a  linear  least-squares  estimator  to  form  an  estimate  of 
the  actual  pixel  value  of  the  true  signal  at  that  pixel  location. 

As  an  alternative,  there  are  other  approaches  which  could  be  used  for 
thresholding  shadows  in  images  with  histograms  that  do  not  have  a  clearly 
bimodal  shape.  One  type  of  approach  involves  iterative  relautation  methods 
[Rosenfeld,81] .  Such  a  method  assigns  probabilities  to  a  pixel's  classification 
and  uses  the  probabilities  of  neighboring  pixels  for  iterative  reinforcement. 
This  iteration  eventually  converges  with  high  probabilities  accumulating  for 
each  class.  It  is  claimed  that  the  method  is  robust.  Such  methods  can  also  be 
extended  to  include  geometric  properties,  not  just  pixel  intensity. 

Another  type  of  approach  is  to  formulate  analytic  criteria  that  the  desirable 
thresholding  is  to  achieve.  The  approach  used  in  [Dunn, 84]  confutes  that 
threshold  which  equalizes  the  probabilities  of  raisclassifying  pixels  in  the 
background  and  object.  This  method  can  also  be  extended,  although  not  as 
efficiently,  to  conpute  a  solution  for  the  minimum  probability  of 
misclassification  error.  The  method  seems  quite  effective  for  highly  speckled 
cind  noisy  imagery  because  of  some  of  the  statistical  independence  assumptions 
made  for  neighboring  pixels  in  (2x2)  windows.  However,  the  analysis  explicitly 
ignores  statistical  contributions  of  pixel  neighborhoods  on  borders  between 
object  and  background.  Therefore,  this  method  may  not  perform  well  on  images 
vdiere  the  background-object  borders  represent  a  significant  fraction  of  the 
total  pixel  content  of  cin  image. 

The  method  of  [Pun, 81]  develops  an  analytic  measure  of  asymmetry  of  the  image 
histogram  based  on  the  information  theoretic  definition  of  entropy.  This  method 
claims  to  have  obtained  some  good  results  on  images  with  quickly  decaying 
histogreims,  such  as  high-pass  filtered  imagery. 

The  approach  used  in  [Perez, 87]  uses  the  illumination-reflectance  model  of 
[Stockham,72]  for  imagery.  In  this  view,  image  intensity  is  the  product  of  a 
reflectance  component  associated  with  the  scene  content  and  an  illumination 
component  which  may  be  spatially  varying. 
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This  method  develops  a  spatially-varying,  adaptive  threshold  function  which 
can  be  implemented  in  raster  fashion.  This  threshold  function  is  approximated 
by  a  highly  differentiable  function  whose  Taylor  series  coefficients  are 
m^eled  by  certain  exponential  functior  3.  This  method  is  intended  for  imagery 
subject  to  slowly  varying  illumination.  Therefore,  this  method  is  not 
recommended  for  use  with  heavily  speckled  SAR  images,  but  may  be  appropriate 
for  passive  optical  imagery. 

For  the  actual  SAR  images  presently  used  in  this  shadow  matching  study,  no 
additional  filtering  was  required  prior  to  examining  the  histogram.  However,  an 
experiment  was  performed  using  the  STAR-2  SAR  image  in  the  upper  left  quadrant 
in  Fig.  4.5  to  condition  the  image  for  feature  extraction.  The  results  of 
filtering  this  image  in  the  upper  left  quadrant  using  the  sigma  filter  with 
values  of  sigma  -  0,  15,  30,  45  are  shown  in  the  other  three  quadreints. 

3. 2. 1.2  Vectorizing  Boundaries 

Vectorization  of  a  feature's  boundary  is  the  process  of  representing  that 
boundary  as  a  connected,  ordered  sequence  of  pixels.  The  vectorization 
procedure  used  in  this  effort  consists  of  three  modules,  and  is  appropriate  for 
binary  images. 

The  first  module  is  the  executive  module.  This  module  scans  an  image  window  for 
binary  region  boundary  pixels  and  determines  whether  the  boundary  region 
intersects  the  image  window  or  is  entirely  within  the  window.  In  either  of 
these  two  mutually  exclusive  cases,  the  executive  calls  either  module  two  or 
three  as  appropriate. 

The  executive  module  scans  the  window  as  follows.  First,  it  examines  each  of 
the  four  borders  of  the  image  window  and  searches  for  boundaries  of  binary 
regions.  Upon  detection  of  such  a  boundary,  the  second  module  is  activated.  The 
second  module  then  tracks  the  boundary  until  it  again  intersects  the  window 
border.  This  module  records  consecutive  boundary  pixel  locations  in  a  special 
data  structure,  as  well  as  in  a  raster  format  array  to  keep  track  of  visited 
locations. 

Next,  the  executive  module  examines  the  interior  of  the  window  line  by  line, 
and  calls  the  third  module  whenever  a  region  boundary  is  encountered  that  is 
not  already  identified  as  a  previously  visited  location  in  the  raster  format 
array.  The  third  module  follows  the  edge  of  a  binary  region  until  it  returns  to 
the  starting  pixel.  This  module,  just  as  the  second  module,  also  records 
consecutive  boundary  pixel  locations  in  a  special  data  structure  and  a  raster 
format  array. 

3.2.2  Matching  Techniques  Using  Shadows  to  Estimate  Sensor  Position 

In  section  3. 2. 2.1,  there  is  a  description  of  some  of  the  algorithms  used  to 
estimate  the  sensor  position  using  measures  of  similarity  between  real  and 
simulated  radar  shadows.  In  section  3. 2. 2. 2,  there  are  some  preliminary  results 
using  some  of  these  methods. 
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3. 2. 2.1  Discussion  of  Algorithms 


Once  a  shadow  region  has  been  thresholded  in  a  real  image,  the  outer  boundary 
must  be  extracted  and  vectorized.  Some  low-pass  filtering  may  be  necessary  if 
there  is  a  great  disparity  in  the  resolutions  of  the  real  imagery  auid  terrain 
model. 

Matching  of  the  shadow  boundaries  in  the  real  and  synthetic  images  can  be 
accomplished  in  a  number  of  ways.  One  type  of  metric  would  seek  to  match  entire 
shadow  regions  in  the  real  aind  synthetic  images.  Such  metrics  are  briefly 
mentioned  in  section  1.3. 3.1,  eind  often  involve  correlation  or  sum  of  aibsolute 
differences.  Because  the  shadow  areas  involve  more  pixels  them  the  shadow 
boundaries,  such  methods  should  be  quite  robust  to  smaller  errors  in  the  data, 
but  less  sensitive  to  small  cheinges  in  the  platform  position. 

Other  methods  for  matching  could  involve  boundary  matching  techniques  (see 
section  1.3. 3.2).  These  methods  would  seem  to  have  the  potential  for  higher 
sensitivity  to  both  platform  position  and  errors  in  the  data  than  area-based 
methods . 

Both  of  these  methods  are  used  for  estimation  of  platform  positioning.  The 
area-based  techniques  are  used  for  acquisition  and  the  contour-based  techniques 
are  used  for  finer  adjustments. 

There  follows  a  discussion  of  the  current  method. 

The  actual  imaging  scenario  is  shown  in  Figure  3.2.  The  real  slant  range  image 
is  created  by  measuring  slant  ranges  r^^j  and  r^^  obtained  from  the  actual 
sensor  position  p  to  the  points  Pj  and  P, .  Then  the  ground  range  image  is 
created  by  projecting  these  slant  ranges  into  the  ground  as  if  the  sensor  were 
at  the  assumed  position  f>.  All  ground  range  measurements  are  relative  to  . 

The  assumption  is  that  p  is  at  altitude  h^^  and  offset  g^^  the  1st  ground  range 
line.  Moreover,  g.^  is  computed  so  as  to  have  slant  range  equal  to  the 
measured  nearest  slant  range  r^^ . 

Here,  the  origin  is  at  the  perpendicular  projection  to  the  ground  of  the 
assumed  sensor  position.  This  is  not  the  origin  of  the  DEM  model.  However, 
all  the  units  can  be  converted  to  DEM  units. 

The  pixel  values  for  Pj  and  p^  in  real  ground  range  image  are: 


"b  - 

-  gAo  )/^gr 

m)  = 

(gA-i 

-  )/^gr 

where 

: 

Ag^.  =  ground  range  resolution 

The  scenario  for  the  simulated  image  is  shown  in  Fig.  3.3.  The  simulated  slant 
range  image  is  created  by  measuring  the  slant  ranges  from  the  simulated  sensor 
position  p^  to  the  model  points  p,  and  p^  .  Denote  these  slant  ranges  by  r^ , 
and  r.  ,  . 

D  1 
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The  simulated  ground  range  image  is  created  by  projecting  the  slant  ranges  r^  3 
and  r j ^  into  the  DEM  ground  plane  from  the  (incorrect)  sensor  position  p, . 

Now  p*  is  at  altitude  \  and  offset  -g^p  from  the  origin  of  the  DEM. 

Denote  these  simulated  ground  reinge  projections  of  Pj  eind  p^  by  g,  3  and  g,  ^ . 

The  pixel  values  for  P3  and  p^  are:  -  (9,37^94)/  where  Ag^  is 
ground  reuige  resolution. 

(1)  An  initial  guess  to  use  for  a  straight  line  approximation  of  the  actual 
flight  path  of  the  sensor  relative  to  the  scene  is  the  assumed  altitude  and 
assumed  ground  range  to  the  first  slant  range  line.  This  first  guess  is  then 
refined  using  shadow  feature  correspondences  between  the  real  image  eind  a 
sequence  of  simulated  images  created  using  refined  position  estimates. 

(2)  Extract  the  shadowed  regions  in  the  real  image  using  the  histogram 
technique  described  in  section  3.2.1,  amd  then  vectorize  the  shadow  boundaries 
as  described  in  section  3.2.2. 

(3)  To  find  correspondingly  shadowed  regions,  use  a  hierarchical  form  of 
normalized  correlation  (see  section  1.3. 3.1)  which  employs  a  sequence  of 
matches  on  a  pyramid  of  reduced  resolution  versions  of  the  real  and  synthetic 
image  taken  along  some  assumed  approximate  flight  path.  Coarse  registration  is 
performed  with  the  coarsest  resolution,  leading  to  finer  adjustments  on  higher 
resolution  versions  of  the  imagery  (see  [Rosenfeld,84 ) ) .  Such  a  scheme  uses  the 
wider  peak  properties  of  the  correlation  metric  to  full  advantage,  since  wider 
peaks  allow  a  wider  "pull-in"  range  for  acquisition  even  when  the  resolution  is 
degraded. 

Such  a  procedure  is  much  more  computationally  efficient  than  exhaustive,  high- 
resolution  correlation. 

(4)  Using  the  vectorized  description  of  shadow  boundaries  in  real  and  synthetic 
images,  form  the  psi-s  representation  [Ballard, 82] .  This  representation  of  a 
cur/e  ideally  gives  tangent  angle  as  a  function  of  arclength  (see  section 
1.3. 3. 2) 

Because  of  the  differences  in  shapes  between  corresponding  real  and  synthetic 
shadows  due  to  terrain  model  errors  and  incorrect  platform  position  estimates, 
the  psi-s  matching  technique  described  in  section  1.3. 3. 2  requires  good  sensor 
position  estimation  in  order  for  the  corresponding  shadows  to  be  similar 
enough . 

Also,  the  psi-s  representation  is  used  to  infer  corresponding  shadow  points  of 
high  curvature  using  the  derivative  of  the  psi-s  representation.  In  this  way, 
this  representation  is  used  as  part  of  an  interest  operator. 

On  the  near  side  of  the  shadow  to  the  sensor,  these  corresponding  points  are 
better  matched  than  on  the  far  side.  This  is  because  on  the  far  side,  the 
shadow's  movement  is  much  more  sensitive  to  errors  in  platform  position  than  on 
the  near  side  where  the  occluding  terrain  elements  lie. 
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At  this  point,  it  is  important  to  recall  that  the  slant  range  to  ground  range 
transformation  of  the  real  SAR  image  was  created  assuming  the  geometry  of  Fig. 
3.2.  This  transformation  was  created  using  the  a  priori  values  of  sensor 
height,  and  range  offset.  Therefore,  whenever  updated  estimates  for  these 
parameters  are  obtained,  they  must  be  used  in  the  creation  of  a  synthetic  sleuit 
rainge  image.  However,  the  a  priori  values  must  be  used  for  the  sleuit  to  grormd 
range  transformation  of  a  synthetic  image.  Otherwise,  a  synthetic  image  could 
be  created  from  exactly  the  correct  sensor  position  and  still  appear  slightly 
different  than  the  real  image. 

(5)  We  now  seek  to  find  that  simulated  sensor  position  p^  that  creates  slant 
ranges  for  Pj  and  p^  which,  when  projected  into  a  ground  ranqe  image  from  p* , 
will  give  the  same  pixel  values  for  Pj  and  p^  as  are  found  in  the  real  ground 


range  image 

( see 

Fig.  3.2  and  3 

.3). 

From  pj  : 

-X3)^ 

+  (h,-h3  )'  =  (g33 

+  gAo)^+h^' 

From  p^  ; 

-^4)^ 

+  (h,-h4  )'=(g54 

+  gAo)'+h^' 

want: 

9,3 

“  9a3~9ao 

9,4 

"  9  A4~9aO 

recall : 

9a  3 

“  m3*^gr+9AO 

9a  4 

-  m/Ag,+g^o 

This  implies 

the 

simultaneous  system: 

-X3  )2+(h3-h3 

=9^  ^  “  Cl 

-Xj  )^+(h5-h^ 

=9^  h;,'  =  C3 

which  is  the  intersection  of  2  circles. 

Another  issue  involves  multiple  measurements  and  the  distribution  of 
measurements.  If  at  least  two  range  measurements  are  available  per  range  line, 
then  the  flight  path  is  determined  for  each  range  line.  This,  of  course, 
ignores  measurement  error.  If  one  could  assume  a  perfectly  straight  and  Ipvel 
flight  path  of  known  direction  then  one  would  only  require  two  lange 
measurements,  not  necessarily  on  the  same  range  line.  Again  measurement  error 
is  ignored. 

A  more  realistic  scenario  would  be  to  assrrme  that  the  direction  of  the  flight 
path  is  known  to  sufficient  accuracy,  that  the  INS  positional  errors  for  the 
imaging  period  give  good  relative  positional  estimates,  and  that  what  is 
required  is  an  offset  estimate  in  crossrange  and  altitude. 
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To  show  how  to  combine  INS  measurements  and  image-derived  range  measurements 
under  these  assumptions,  we  first  demonstrate  how  to  use  multipie  measurements 
for  estimating  a  constauit  horizontal  flight  path  of  known  orientation.  This 
technique  borrows  from  a  method  used  in  sonar  [Smith,  87]  for  locating  a  source 
given  an  over-determined  system  of  simultcuieous  rauige  measurements.  This  method 
will  then  be  generalized  for  the  use  of  INS  measurements. 

Let: 

Xg  -  sensor  position 

-  known  DEM  locations,  i-1  _ _ N 


Using  the  N  geographic  DEM  points,  one  can  generate  ein  over-determined  set  of 
nonlinear  equations.  However,  using  range  differences  emd  one  known  range 
value,  one  caui  generate  ein  over-determined  set  of  linear  equations. 

From  the  definition  of  c^^  we  have: 

|2-2x%Xi  +  |X^  12-0^2 

Now  let: 

1 

-  l^i  I 

Using  these  definitions  and  that  for  and  letting  one 

of  the  Xj^  be  the  origin,  we  have  for  i-2,  ...,  N 

ir:,+ 

This  implies: 

0  -  r^^  -  ACi^-2r,  ACii-2x'^^x7  i  -  2,  ...,  N 

Note  that  the  unknown  sensor  position  x^  now  enters  linearly  into  this  system  of 
equations,  even  though  range  is  a 

nonlinear  function  of  x^  . 

Let  =  "equation  error"  for  i*^^  equation,  then: 

-  Ac^j-2rg  ACj^^-2j?'^X5  i  =  2,  ...,  N 
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Using  vector  eund  matrix  notation,  this  system  of  equations  cam  now  be 
summarized: 

T  -  V  -  2rj  ^  -  2MXj 

where: 

£  “  (  Ej  /  •  •  •  Eu 
V  -  (  Tj  2  -  ‘^^2  1  '  •  •  •  ' 

M  -  [Xj  Zj" 

•Xn  Zn- 

Ac  "  (  AC2  ^  ,  •  •  •  ,  1 

This  over-determined  system  cam  be  solved  to  yield  the  minimum  least  squares 
value  for 

by  the  use  of  the  generalized  inverse: 

X,  -  1/2  MT  (v  -  2r,  Ac) 

Such  an  expression  can  even  be  generalized  to  include  a  weighting  matrix  W  for 
the  measurements: 

X,  -  1/2  (M^WM)-^  MFW  (v-2r,Ac) 

Now,  we  generalize  this  result  to  incorporate  the  use  of  INS  measurements  and  we 
again  assume  that  what  is  required  is  an  offset  estimate  in  cross  range  and 
altitude. 

Let; 

{t.  }  be  the  time  points  corresponding  to  a  range  measurement  with 
respect  to  an  identified  DEM  position. 

{pCt^)}  be  the  set  of  INS  measurements  for  the  time  points  {t^}  i  » 
1,  ...,  N 

Xg  be  the  desired  INS  offset 
Then  make  the  modified  definition: 

Ci  -  lp(ti  )  -  )  I 

but  c^  =■  |x^  -  (x^  -  p(ti  ) )  I 
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Therefore,  one  can  use  the  previous  generalized  inverse  solution  to  obtain  the 
INS  offset  but  with  the 

replaced  by  (x^  -  p(t^ )) 

Of  course,  amother  smaller  error  will  occur  because  the  association  between 
image  and  model  points  is  dependent  on  occluding  points  for  shadows.  If  the 
sensor  position  cheuiges,  the  position  of  occluding  points  chauiges  slightly. 
However,  this  is  a  second-order  effect,  and  will  be  ignored  because  the  process 
will  be  repeated  by  iteration  for  the  new  estimated  sensor  position  P'  . 

(6)  Using  the  newer  approximation  for  the  flight  path,  iteration  can  now  be 
performed  on  the  previous  steps.  The  convergence  should  be  rapid. 

Of  course,  the  problem  of  resection  for  any  sensor  is  ill-conditioned  if  the 
control  points  are  not  well-distributed.  For  the  case  of  SAR,  some  exeunples  are 
given  in  Fig.  3.42  -  3.45  for  the  variation  in  the  location  of  the  intersection 
of  two  reuige  circles  as  a  function  of  perturbations  of  the  center  separations. 

Another  relevant  problem  is  the  degree  of  precision  required  for  the  DEM  to 
achieve  a  given  accuracy  for  shadow  matching.  Some  numerical  results  on  this 
issue  are  shown  in  Fig.  3.46  -  3.49. 

3. 2. 2. 2  Results 

In  order  to  test  some  of  the  concepts  in  the  previous  section,  one  test  that  was 
performed  was  to  vary  the  assumed  altitude  of  the  flight  path  (assuming  that  the 
flight  path  is  a  straight,  horizontal  line)  and  create  synthetic  images 
corresponding  to  these  offsets.  The  real  SAR  image  is  in  Fig.  3.4.  Synthetic 
images  for  a  number  of  different  altitudes,  but  correct  ground  range  offsets, 
are  shown  in  Fig.  3.11  -  3.17. 

The  use  of  normalized  correlation  metric  for  acquisition  matching  of  shadowed 
regions  is  shown  in  Fig.  3.18  -  3.24.  A  large  window  was  picked  out  in  the  real 
image  which  contained  a  shadow.  This  was  established  using  the  thresholding  in 
step  (1)  of  the  algorithm.  Because  of  an  initial  positional  approximation 
available  from  the  INS,  it  is  possible  to  select  a  small  region  within  the 
simulated  images  which  will  be  known  to  be  within  this  larger  region  in  the  real 
image.  Then,  the  normalized  correlation  metric  was  used  to  find  that  patch 
within  the  real  image  which  corresponded  to  each  of  these  patches  in  the 
synthetic  images  corresponding  to  the  different  flight  altitudes.  The  best  match 
was  achieved  using  a  hypothesized  value  of  10.0  km. 

The  use  of  the  sum  of  differences  metric  for  matching  shadowed  regions  is  shown 
in  Fig.  3.25  -  3.29.  The  same  procedure  was  employed  as  for  the  images  in  Fig. 
3.18  -  3.24. 

The  vectorized  shadow  boundaries  corresponding  to  different  flight  altitudes  are 
shown  in  Fig.  3.31  -  3.37. 

The  computed  flight  path  altitude,  for  which  the  best  match  using  the  normalized 
correlation  metric  was  achieved,  was  10.0  km.  This  agrees  with  the  flight  log. 
Fig.  3.41  shows  the  values  of  the  metrics  used,  i.e.  normalized  correlation  and 
sum  of  differences,  on  both  intensity  and  binarized  images. 
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Fig.  3.38  shows  a  windowed  area  containing  shadows  used  for  the  correlations, 
along  with  the  points  identified  from  their  high  curvature  in  their  psi-s 
representation.  The  same  is  shown  in  the  corresponding  synthetic  image  shown  in 
Fig.  3-39.  Finally,  Fig.  3-40  shows  the  correspondences  formed  between  these 
identifiable  points  in  the  real  and  synthetic  images.  These  corresponding  points 
are  used  to  implement  steps  (5)  and  (6)  of  the  algorithm  described  in  the 
previous  section. 

The  particular  SAR  image  shown  in  Fig.  3.4  has  its  radar  shadows  clustered 
together  in  the  range  direction  since  the  single  mountain  range  in  the  field  of 
view  is  roughly  parallel  to  the  flight  path.  Therefore,  using  shadows  alone  in 
this  case  will  lead  to  ill-conditioned  resection  unless  other  control 
information  is  used  that  is  more  spread  out  in  the  range  direction.  In  general, 
various  combinations  of  control  information  from  both  terrain  auid  features  will 
be  required  to  form  well-conditioned  resections. 

3.2.3  Analysis  of  Robustness 

Using  the  previous  imagery  and  the  positional  acquisition  technique  described 
aibove,  several  initial  estimates  for  offset  were  attempted.  The  limits  on 
initial  errors  are  discussed  below. 

In  one  case,  the  correct  sensor  position  was  about  10  km.  in  altitude,  and  14 
km.  in  reinge  to  near  range  offset.  The  correlation/matching  technique  failed  at 
initial  estimates  for  altitude  in  excess  of  12  km.  and  lower  thcui  7.5  km.,  ie. 
for  about  a  25%  relative  error. 

In  another  case,  where  the  correct  ground  range  offset  was  9.76  km.,  all  initial 
estimates  within  5  km.  to  20  km.  converged  to  the  correct  solution.  These  two 
upper  and  lower  values  represent  only  the  limits  of  testing. 

Azimuthal  offsets  of  up  to  27  km.  were  tested,  and  convergence  occured  in  all 
cases. 

Clearly,  the  above  comments  only  represent  empirical  results  using  a  small  data 
set.  A  more  complete  analytical  analysis  must  be  performed  to  better  understand 
the  robustness  of  these  procedures. 

3.3  Use  of  Terrain-Induced  Shading  for  Control 

The  use  of  shading  information  in  a  real  and  synthetic  optical  image  for 
registration  of  image  and  terrain  model  has  been  done  previously  [Horn, 78]. 
However,  the  transformation  tHat  was  sought  was  not  a  perspective,  but  just  a 
global  rotation  and  translatioti  since  the  imaging  scenario  was  orthographic 
projection. 

The  use  of  SAR  imagery  presents  seme  differences  from  the  case  reported  in 
[Horn, 78).  SAR  imagery  has  more  artifacts  such  as  noise  and  speckle  than  optical 
imagery,  and  radar  backscatter  has  a  more  directional  dependency.  Also,  SAR 
imagery  is  side-looking  rather  than  down-looking. 
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However,  the  above  approach  can  still  be  used  for  selected  regions  in  the  SAR 
case  with  some  modification.  Step  (2)  of  the  algorithm  in  the  previous  section 
3.2.2  used  region-based  correlation  to  approximately  identify  corresponding 
shadowed  regions  previously  extracted  in  step  (1).  By  using  the  model,  an 
algorithm  could  select  other  favorable  terrain  for  matching.  This  is  an  example 
of  "model-based"  matching. 

For  excimple,  fairly  level  regions  which  are  near  the  tops  of  elevated  terrain 
would  not  cause  rauige  compression,  and  may  have  good  contrast  with  possible 
shadows  on  the  other  side  of  the  crest. 

Some  local  histogram  re-mapping  may  also  be  useful  when  matching  cemdidate 
regions  to  reduce  the  effect  of  variations  in  scattering  cross  sections. 

This  approach  may  also  be  investigated  in  the  siibsequent  Options. 


45 


Sensor 
I  Pusiiion 


1> 


,s 

c/} 

r3  w  c 
-■&-2 
rn  •«  5 

-o'  =i3 

2S 

O 

u. 


46 


Figure  3.1(b) 
Vectors  corresponding  to 
the  surface  gradient,  and 
the  surface  normal. 
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Vectorized  shadow  boundries  in  simulated  image 
Brazeau  Range  Area,  Alberta,  Canada 
Sensor  Height;  12kra 
Range:  15.49km 

Maximum  Height  in  model:  2.5km 
DEM  Resolution:  60m  x  60m  x  Im 


Figure  3.36 


Vectorized  Shadow  Boundaries  in  actual  SAR  image 
Brazeau  Range  Area,  Alberta,  Canada 
Approximate  Altitude:  10km 

Approximate  Range:  14km 
SAR  imagery  by  INTERA  Technologies 


Figure  3.37 


Matched  Points  in  Brazeau  Range  Simulated  Radar  Image 
Flight  Altitude:  10km 
Correlation  Window:  287,313,50,50 


Figure  3.39 
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Figure  3.42 


Intersections  of  Two  Circles  with 
Horizontal  Perturbation  to  Far  Point 
Center  Separation  Approx.  0,1R 
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Figure  3.43 
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Figure  3.45 
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4.0  Generation  of  Control  Infonnation  Using  Feature  Data 

Some  of  the  methods  described  in  section  1.3. 3. 2  on  bovindary  matching  are 
relevant  to  matching  the  outlines  and  shapes  corresponding  to  feature  data. 
Region  matching  methods  may  also  be  applicable,  but  may  suffer  some  performance 
problems  if  intensity  distortions  are  not  predictable.  Of  course,  the 
descriptions  of  feature  data  will  be  the  coordinates  of  the  endpoints  of 
polygonal  approximations. 

The  appearance  in  an  optical  image  of  such  a  polygonal  approximation  will  be 
modified  by  a  perspective  transformation.  A  ground  range  SAR  image  will  have 
some  distortions  if  such  features  lie  on  terrain  relief. 

However,  if  these  features  are  not  distorted  by  terrain,  then  the  matching 
techniques  described  in  section  1.3. 3. 2  are  applicable.  Because  flight  path 
orientation  is  expected  to  be  known  within  a  degree,  only  translation  must  be 
estimated  for  a  ground  range  SAR  image. 

Perspective  must  be  estimated  and  accounted  for  in  an  optical  image.  Of  course, 
the  perspective  distortion  for  an  object  at  a  hypothesized  distsmce  from  the 
optical  sensor  is  predictable  from  this  hypothesized  offset  distauice.  Knowing 
the  approximate  height  of  the  sensor  platform,  the  sensor  orientation,  and  a 
region  of  pixels  that  is  being  searched  for  a  given  feature,  one  can  approximate 
the  distance  of  the  feature  from  the  sensor.  Then  one  can  approximately  predict 
the  perspective  distortion.  This  distortion  should  not  vary  appreciably  over  the 
field  of  view. 

Searching  for  a  particularly  distinctive  feature  in  the  feature  database  is  a 
plausible  strategy,  as  is  discussed  in  section  7  concerning  rules  for  choice  of 
features.  However,  another  strategy  is  to  extract  features  from  the  imagery  and 
try  to  very  efficiently  find  the  corresponding  features  in  the  database.  This 
latter  approach  will  be  discussed  for  future  work. 

The  procedure  that  was  initially  explored  for  feature  matching  involved 
searching  efficiently  for  one  particular  feature.  The  steps  involved  in  using 
the  NHAP  photography  were  pre-processing  the  image  to  reduce  the  amount  of 
micro-edges  due  to  noise  and  micro-structure,  ^md  using  a  partial  matching 
scheme  to  reduce  the  computational  conplexity  of  examining  trial  positions  of 
the  feature  (see  Fig.  4.16). 

A  more  detailed  summary  of  the  steps  follows: 


1. 

^ply  Sobel  edge  operator  [Pratt, 78)  to  image. 

2. 

Threshold  result  at  +2  sigma  distributional  gray 
to  retain  major  edges. 

value 

(see  Fig. 

4.17) 

3. 

Delete  isolated  pixels. 

4. 

Convert  selective  feature  point  coordinates 
coordinates . 

from 

UTM  to 

image 

5. 

Subsample  the  rasterized  feature  (see  Fig.  4.16). 
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6.  Correlate  subsampled  feature  pixels  with  pre-processed  image  derived  at 
step  (3),  skipping  offsets  where  a  first  pixel  of  the  subset  does  not 
overlay  a  bright  pixel. 

7.  Retain  locations  of  top  three  correlation  maxima,  and  test  the  full 
feature  pattern  at  these  points  to  obtain  best  estimate. 

8.  Use  localized  edge  detection  to  obtain  better  estimates  for  tie  points 
in  feature  and  original  image. 

9.  Use  tie  points  to  determine  resection. 

Various  combinations  of  pre-processing  steps  will  be  investigated.  Included 
among  these  will  be  resolution  reduction  to  reduce  amount  of  microstructure, 
edge  detection,  and  thresholding.  Some  of  these  steps  are  shown  in  Fig.  4.2- 
4.15. 

SAR  image  feature  matching  used  STAR-2  imagery  of  the  same  general  Brazeau 
region  as  for  shadow  matching,  described  in  sec.  3,  with  reduction  of  resolution 
to  512  X  512  pixels. 

Again,  various  combinations  of  pre-processing  operations  can  be  used  to  improve 
the  efficiency  of  the  feature  matching  process  for  a  particular  selected 
feature.  Combinations  of  the  following  steps  will  be  explored: 

R  «  resolution  reduction 

L  «  Lee's  sigma  filter  (see  sec. 3. 2. 1.1) 

S  -  Sobel  edge  detection 

T  -  thresholding  using  bimodality  of  processed  image  histogram 
n  (digit)  -  thresholding  at  +n  std.  dev.  gray  value 
V  -  vectorization  of  binary  image  boundaries 
The  following  combinations  of  these  steps  are  expected  to  be  exzunined: 


LSI 

LS2 

LS3 

LTS 

LTSl 

LTS2 

LTS3 

LTV 

RSI 

RS2 

RS3 

RTS 

RTV 

RLSl 

RLS2 

RLS3 

RLTS 

RLTV 

LTVR  LTSR 

Another  approach  that  will  be  explored  involves  major  edge  extraction  from  the 
image  to  obtain  significant  boundaries  that  will  be  used  to  selectively  search 
the  feature  database.  Relaxation  methods  have  been  used  to  some  extent  for 
obtaining  such  major  contours  [Rosenfeld,82] ,  and  will  be  explored  in  the  next 
effort. 


In  this  approach,  the  feature  database  is  only  selectively  searched  becaiise  the 
dat2Q>ase  has  been  orgcuiized  in  geometric  feature  groupings  offline.  Sucli 
geometric  groupings  are  eunalogous  to  hash  coding  addresses  in  ordinary 
databases,  and  thus  the  term  "geometric  hashing"  [ Schwarz , 86 ] ,  [Kalvin,87]  (see 
sec.  1.3.3).  Another  unexplored  representation  uses  fractuals  for  feature 
representation  I  Barnsley,  88  ] .  This  representation  has  been  used  for  image 
compression,  but  not  for  model-based  recognition.  Such  approaches  for  reducing 
database  search  will  be  investigated  in  the  follow-on  Options. 

The  approach  just  described  is  distinct  from  the  previous  approach  which  first 
selects  a  feature  from  the  database  and  then  searches  the  image  for  it. 
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Brazeau  Range  Aerial  Image,  Alberta,  Canada 
NHAP  Photography  Digitized  at  USC 
Frame  176 


Figure  4.1 
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Figure  4.5 
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star-2  Image  of  Brazeau  Range  Area 
Preprocessing;  Lee  sigma  filter 
Sobel  edge  image 


Threshold  at  mean  +  sigma 


Figure  4.7 


STAR-2  Image  of  Brazeau  Range  Area 
Preprocessing:  Lee  sigma  filter 

Sobel  edge  image 
Threshold  at  mean  +  3  *  sigma 

Figure  4.9 
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STAR-2  Image  of  Brazeau  Range  Area 
Preprocessing:  Lee  sigma  filter 

Threshold  by  bimodal  histogram 
Sobel  edge  image 

Figure  4.10 
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STAR-2  Image  of  Brazeau  Range  Area 
Preprocessing;  Lee  sigma  filter 

Threshold  by  bimodal  histogram 
Sobel  edge  image 
Threshold  at  mean  +  sigma 

Figure  4.11 


of  Brazeau  Range  Area 
:  Lee  sigma  filter 

Threshold  by  biroodal  histogram 
Vector i zation  of  edges 


3TAR-2  Image  of  Brazeau  Range  Area 
Preprocessing:  Lee  sigma  filter 

Threshold  by  bimodal  histogram 
Vec to r i za t i on  of  edges 
2X  resolution  reduction 


Fioure  4.14 
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Figure  ^.16 


5.0  Generation  of  Control  lnfonnatic»i  Using  Image  Segment  Data 

Image  segments  or  "chips"  are  another  way  to  describe  a  map  feature.  Such  a 
segment  is  an  image  of  the  given  feature  as  well  as  some  of  the  surroiinding 
area.  One  way  to  find  the  given  feature  in  the  sensed  image  is  to  attempt  to 
register  the  segment  with  some  region  of  the  sensed  image. 

Another  way  would  be  to  also  have  a  description  of  the  given  feature, such  as  a 
polygonal  approximation,  along  with  the  image  segment.  Then  a  combined  boundary 
match  and  regional  match  approach  would  be  possible. 

VEXCEL  has  created  a  particularly  efficient  algorithm  to  estimate  relative 
translations  of  a  pair  of  opposite-side  satellite  SAR  images.  VEXCEL  has  also 
developed  a  boiondary-based  algorithm  for  registering  objects  which  have 
undergone  translation  and  rotation  between  two  images.  These  two  algorithms  are 
described  in  section  1.3. 3. 2. 

Another  approach  for  estimating  relative  translation  of  one  image  relative  to 
another  is  similar  to  phase  correlation,  but  is  possibly  more  efficient. 
Briefly,  this  proposed  method  is  as  follows: 

(1)  Perform  some  type  of  high-pass  filtering  of  both  images  to  enhance  the 
edges  and  de-emphasize  the  slowly  varying  regional  content. 

(2)  Conpute  the  FFT  of  both  images. 

(3)  For  each  conplex  frequency  f^^’  in  image  1  and  the  corresponding  complex 
frequency  fi^>  in  image  2,  form  the  normalized  difference 

^ M/(2nifj'* ’ )  where  are  the  corresponding  phases  (formed  from 

tie  quotients  of  the  magnitudes  of  the  imaginary  and  real  parts  of  the 
complex  amplitude  of  a  frequency). 

(4)  Analogous  to  the  Hough  transform,  form  clusters  of  the  values  of  these 
normalized  phase  differences.  By  the  Fourier  Translation  [Champeney,  73] 
Hieorem,  a  relative  phase  difference  for  the  same  frequency  in  two  images 
corresponds  to  a  relative  translation  of  this  frequency  in  the  two  images. 
Any  trauislation  of  a  substantial  anount  of  feature  data  will  lead  to  a 
substantial  number  of  corresponding  frequencies  shifted  by  a  common  phase 
difference  ^unount.  This  common  phase  difference  will  be  observable  as  a 
cluster  of  values  in  the  list  of  corresponding  phase  differences. 

This  method  has  potential  for  sub-pixel  registration,  since  the  phases  arise  as 
quotients  of  the  magnitudes  of  imaginary  and  real  frequency  amplitudes.  These 
quotients  are  not  limited  to  integral  values. 

This  method,  along  with  other  image-image  matching  schemes,  will  be  investigated 
during  the  Options  to  the  Base  Contract. 


Ill 


6.0  Generation  of  Control  Infonpaticm  Using  Satellite  Epheneris  Data 

Registration  of  orbital  SAR  imagery  to  maps  is  difficult  because  of  range 
distortions  due  to  terrain  relief  and  azimuthal  distortions  due  to  earth 
rotation  and  pointing  errors.  In  the  past,  pointing  errors  have  involved  a 
linear  skew  for  imagery  such  as  SEASAT  and  SIR-B.  This  skew  was  removable  by  the 
use  of  control  points  and  rubber  sheet  interpolation.  However,  these  pointing 
errors  will  be  a  problem  for  SIR-C  aiui  X-SAR  imagery  because  of  the  nonlinear 
distortions  associated  with  their  associated  larger  doppler  shifts 
(Curlander,87] . 

Other  errors  and  their  sources  for  satellite  SAR  are  [Curl2mder,82} : 
o  Rcinge  errors 

geoid  errors 

undetected  delays  in  data  collection  hardware 
gro\ind  range  nonlinearities 

incidence  angle  varies  across  swath,  so  ground  resolution  of 
each  time  sample  is  not  uniform 

incidence  angle  is  function  of  terrain 

o  Azimuth  errors 

synchronization  errors  between  ephemeris  data  and  returning  radar 
echoes 

GSFC  tapes  have  ephemeris  data  for  1  m  intervals 
SEASAT  SAR  processing  interpolates  these  for  1  ms  intervals 
missing  azimuthal  lines 

azimuth  skew  because  doppler  shift  is  function  of  latitude  and 
position  within  swath 

Azimuth  skew  can  be  dealt  with  by  more  complex  processing  with  respect  to  zero 
doppler,  or  by  post-processing  deskewing  using  control  points. 

Also,  the  look  angle  tends  to  amplify  ground  location  errors,  since; 


AEj  -  AR/sina 
where; 

a  «  radar  angle  of  incidence  from  vertical 
AR  ■  range  offset 

AEj  -  location  error  on  ground  in  range  direction 
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Therefore,  reducing  this  look  angle  will  reduce  this  ground  range  error  effect. 

Some  of  the  refinements  to  the  process  of  associating  an  earth  location  with 
pixels  are  used  in  a  process  called  "geocoding"  [Curlander,87] : 

o  More  accurate  geoid. 

o  More  accurate  accounting  of  delays  in  timing  calibration  of  r2mge 

measurements. 

o  Taking  into  account  the  SAR  processing-induced  skewing  of  iso-range 

lines  to  more  correctly  conpjte  slant  range  for  a  given  pixel. 

o  Use  of  doppler  function  that  takes  into  account  earth  rotation 

o  Use  of  doppler  plane  at  center  of  radar  beam  as  a  function  of  time 
instead  of  pointing  direction  of  euitenna. 

o  Use  of  terrain  model. 

More  accurate  analysis  of  the  doppler  function  for  the  satellite  case  is 
described  in  (Li, 85],  where  it  is  shown  that  there  are  two  dominemt  terms  for 
the  matched  filter  doppler  function:  doppler  centroid  and  doppler  frequency 
rate. 

The  doppler  centroid  term  is; 

f^c  -  (-2,  .  R,  vJ/(  >Si) 

Errors  in  performing  the  SAR  correlation  with  errors  in  this  terra  of  the  matched 
filter  will  lead  to  degradations  in  the  signal  to  noise  ratio  and  the  signal  to 
ambiguity  ratio. 

The  doppler  frequency  rate  term  is: 

fDR  -  2(v,  •  V,  -  R,  .A,,  )/(  XR) 

Errors  in  this  terra  for  correlating  the  raw  SAR  data  lead  to  de-focusing. 

Here: 

^  (t)  -  position  of  target  on  Earth  surface  at  time  t 
(t)  ■  position  of  spacecraft  radar  at  time  t 
R,  -  It  (0)  -  R.,  (0) 

V.  (t)  -  velocity  of  target  on  Earth  surface  at  time  t 


Vgg.  (t)  -  velocity  of  spacecraft  radar  at  time  t 

V,  -  Vt  (0)  -  V.,  (0) 

-  acceleration  of  spacecraft  radar  at  time  t 
X  -  radar  wavelength 

Higher  quality  attitude  and  ephemeris  information  can  lead  to  improvements  in 
processing  the  raw  SAR  data  using  these  frequency  terms.  However,  methods  are 
also  discussed  in  [Li, 85]  for  estimating  these  two  terms  from  the  raw  radar 
return  data. 


7.0  Rule-Based  Issues  for  Generation  of  Control  information 
Five  sources  of  control  information  are  described  in  this  report: 
o  INS 

o  GPS 

o  Terrain-induced  shadows 

o  Terrain-induced  shading 

o  Pleunimetric  features 

What  is  required  is  an  overall,  top-down,  coherent  strategy  for  using  these 
various  sources  and  their  combinations  under  the  appropriate  circtmnstances.  All 
of  these  sources  have  advantages  and  disadvantages,  requirements,  and  various 
degrees  of  achievable  control  accuracy. 

One  way  to  formulate  the  expertise  required  to  correctly  exploit  these  multiple 
sources  is  by  the  use  of  a  rule-based  system.  Other  options  for  incorporating 
heuristic  knowledge  include  the  use  of  frame-based  methods  [Minsky, 75],  which 
use  a  template  approach  for  representing  "typical"  situations. 

There  are  some  potential  problems  when  dealing  with  rules.  A  rule-based  approach 
can  be  difficult  to  inclement  for  situations  involving  complex  relationships 
(Mettrey,87} .  Such  systems  could  potentially  become  difficult  to  understand  when 
there  is  a  large  number  of  rules,  unless  the  rule  set  is  hierarchically 
organized. 

However,  rules  do  provide  a  format  for  incorporating  knowledge  that  is 
relatively  sinple  eund  easy  to  use.  Moreover,  the  approach  that  will  be  followed 
will  be  to  use  a  top-down  approach  for  creating  and  organizing  the  rule  set. 
Therefore,  the  rule  set  shoudd  be  understemded^le  even  as  it  grows. 

Such  a  system  would  be  somewhat  misnaned  as  an  "expert  system  ",  since  an 
argument  may  be  given  that  vhenever  one  is  operating  strictly  according  to 
niles,  one  is  really  operating  at  the  level  of  the  advanced  beginner.  A  real 
expert,  on  the  other  hand,  is  able  to  create  new  rules  based  on  experience, 
modify  and  find  appropriate  exceptions  to  existing  rules,  euid  can  make  novel 
hypotheses . 

Nevertheless,  the  output  from  such  a  rule-based  system  is  potentially  very 
helpful  in  average  situations  which  would  not  require  cuiy  extraordinary 
interpretation  by  a  human  expert  in  radargrammetry.  We  shall  call  such  a  system 
a  "nile-based  assistant". 

Such  aui  rule-based  assistant  would  decide  which  procedural  modules  to  execute 
and  would  evaluate  their  outputs  using  its  knowledge  base  consisting  of  rules. 
These  rules  are  based  on  heuristics  and  "rules  of  thumb"  commonly  enployed  by 
practicing  experts. 
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In  this  sense,  the  value  of  a  rule  should  be  its  applicability  in  a  wide  variety 
of  cases.  The  existence  of  exceptional  cases  should  not  detract  from  the  value 
of  a  particular  rule. 

Rules  with  negative  conclusions  are  important  also.  For  example,  a  rule  which 
cautions  against  picking  a  certain  feature  under  a  particular  circumstance  is  as 
useful  as  one  vdiich  positively  chooses  it  under  different  circumstances.  Such 
reasoning  allows  the  system  to  eliminate  possibilities. 

Ihere  is  often  a  "gray  area"  between  positive  and  negative  rules  v^erein  it  is 
not  clear  that  a  particular  feature  is  suitable  or  not.  Such  ambiguities  result 
because  rules  that  eliminate  a  candidate  feature  are  not  necessarily  converses 
of  rules  that  positively  choose  a  candidate. 

The  result  is  that  most  rules  do  not  give  necessary  and  sufficient  conditions 
for  choosing  a  particular  feature  for  control.  Therefore,  there  is  always  the 
possibility  of  inserting  new  rules  that  are  inconsistent  with  the  existing 
knowledge  base.  This  is  always  a  potential  problem  for  large  rule-based  systems. 

However,  such  occurences  should  be  relatively  rare  if  the  knowledge  base  is 
suitably  organized  in  a  top-down  hierarchy.  Then  many  rules  are  often  partial 
converses  of  other  rules. 

Although  the  present  study  is  not  concerned  with  the  explicit  construction  of 
such  a  rule-based  system,  the  organization  and  content  of  the  knowledge-base 
will  be  discussed.  The  content  of  this  knowledge  is  first  expressed  in  meta¬ 
rules,  Tidiich  will  be  subsequently  translated  into  actual  coded  n^es  in  a  non- 
proce^ral  language.  Some  of  these  roeta- rules  appear  in  section  7.!z. 

This  rule-based  effort  has  just  begun  under  the  last  part  of  the  Base  Contract, 
and  will  continue  during  the  additional  Options. 

7.1  Overview  on  Rules  for  image  Control 

Regardless  of  sensor  type,  the  output  of  a  rule-based  system  for  generation  of 
image  control  must  determine  the  following  (see  Fig.  7.1): 

o  Choice  of  control  type 

o  Computation  of  control 

o  Estimation  of  control  accuracy 

o  Conflict  resolution 

The  categories  of  control  information  include  INS,  GPS,  terrain-induced  shadows, 
terrain-induced  shading,  and  planimetric  features  (see  Fig.  7.2).  A  selection  of 
a  single  type  or  multiple  types  of  control  would  be  made. 

Choices  may  also  be  appropriate  within  a  control  type.  For  example,  if 
planimetry  is  chosen,  then  the  question  of  which  planimetry  features  to  search 
for  first  is  in^xirtant.  These  selections  are  based  on  both  image  formation  and 
image  processing  considerations.  This  module  is  rule-based. 
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Image  Control  must  be  computed  from  the  control  type  or  types  selected.  This 
module  is  entirely  procedural,  euid  represents  con^itations  such  as  in  sec. 
3.2.2. 

An  estimate  is  also  required  to  determine  the  potential  accuracy  of  control 
given  the  distribution  of  control  features.  This  estimate  may  or  may  not 
indicate  the  necessity  for  obtaining  additional  control  information.  This  module 
is  also  procedural. 

Finally,  conflicts  or  ambiguities  may  arise.  Some  capability  for  resolution  of 
such  problems  is  required  for  deciding  to  obtain  additional  control  or  for  re¬ 
evaluating  the  accuracy  of  obtained  control.  This  module,  shown  in  Fig.  7.1, 
would  resemble  a  diagnostic  system  and  is  rule-based. 

Expanding  the  planimetry  box  in  Fig.  7.2  shows  the  organization  of  how 
planimetry  is  chosen  and  evaluated  as  shown  in  Fig.  7.3.  Two  types  of  reasoning 
are  used  for  evaluating  a  potential  feature. 

One  way  is  to  use  "first  principles"  from  the  basic  physics  of  electromagnetics 
^ich  qualitatively  evaluate  the  strength  of  the  signal  response.  Another  way  is 
to  use  "rules  of  thumb"  from  practicing  photo-interpreters  and  image  processors. 

In  the  area  of  image  formation,  such  ei^irical  rules  are  also  ultimately  related 
to  physical  principles.  However,  they  may  be  easier  or  more  convenient  to  state 
in  a  single  heuristic  form,  rather  thaui  as  a  consequence  of  a  combination  of 
physical  principles. 

In  the  area  of  image  processing,  enpirical  considerations  are  important  also. 
Some  features  are  confiutationally  easier  to  locate  than  others.  Rules 
representing  such  choices  would  also  be  best  categorized  as  "rules  of  thumb". 

Finally,  there  exists  a  module  for  forming  composite  evaluations  of  features 
based  on  multiple  individual  evaluations. 

Fig.  7.4  shows  the  organization  of  the  rationale  for  employing  basic  physical 
principles  to  the  evalxiation  of  planimetric  features.  The  relevant  physical 
characteristics  clearly  pertain  either  to  the  object  to  the  sensor  scenatiu. 
The  sensor  scenario,  shown  in  Fig.  7.6,  relates  to  SAR,  optical,  or  IR.  In  Fig. 
7.7,  the  SAR  scenario  concerns  effects  due  to  wavelength,  polarization,  spatial 
resolution,  the  direction  of  the  flight  path  relative  to  the  object,  amd  the 
signal  power.  The  characteristic  properties  of  the  feature  object  important  to 
imaging  are  shown  in  Fig.  7.8,  and  concern  size,  surface  roughness,  and 
dielectric  constant. 

Rules  of  thumb  are  partitioned,  as  in  Fig.  7.5,  into  considerations  involving 
the  imaging  formation  scenario,  as  shown  in  Fig.  7.9,  or  issues  related  to  the 
con^jutational  complexity  of  image  processing  procedures  for  extracting  the 
feature  from  the  image. 

An  observation  that  can  be  made  concerning  such  a  top-down  developnent  of  rules 
is  that  the  rules  become  more  specialized  and  detailed  as  a  function  of  depth  in 
the  structure  chart  represented  by  Fig.  7.1  to  7.9. 
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7.2  Meta-Rules 


To  organize  the  development  of  meta-niles  as  described  above,  an  indexing  scheme 
will  be  helpful  which  immediately  identifies  the  rule  as  belonging  to  a 
particular  hierarchical  level  as  shown  in  Fig.  7.1  to  7.9.  Therefore,  rules 
vAiich  decide  vdiich  category  of  control  to  select,  as  in  Fig.  7.2,  will  be 
labeled  with  a  prefix  ”2.",  followed  by  an  index  nimiber  which  distinguishes 
rules  at  that  level. 

Meta-rules  are  stated  in  the  form:  If  (hypothesis)  then  (conclusion).  The 
notation  "HM1.N2"  refers  to  the  hypotheses  of  Rule  N1.N2. 

The  accuracy  estimation  module,  shown  in  Fig.  7.1,  will  be  denoted  by  "AM".  This 
module  is  a  procedural  function  of  its  inputs,  l^ese  inputs  can  be  measurements 
from  any  combination  of  control  types.  Another  input  would  be  the  required 
accuracy  E.  Given  these  inputs,  the  function  AM  will  then  decide  vdiether  or  not 
the  input  measurements  support  the  required  control  accuracy. 

As  an  exan^le,  we  denote  v^ether  the  estimated  accuracy  using  planimetry, 
shadow,  and  accuracy  inputs  is  acceptable  or  not  by  AM(plamimetry, shadows, E)  - 
Yes  or  -  No. 

Another  procedural  module  with  inputs  is  the  control  computation  module,  denoted 
cr(  ).  The  input  could  be  of  more  than  one  control  type.  For  example,  the 
notation  CT( planimetry, shadow)  indicates  that  control  should  be  confuted  for  an 
image  with  planimetry  amd  shadow  data. 

Finally,  the  notation  GPS  «  Yes  or  •  No  indicates  whether  or  not  GPS  data  is 
available  in  a  given  scenario.  However,  in  order  to  make  the  notation  tractable, 
an  abuse  of  notation  will  be  permitted.  Terms  such  as  "planimetry"  will  be  used 
to  indicate  an  input  type  for  a  functional  module  such  as  AM  or  CT,  as  well  as 
to  indicate  a  logical  variable  (which  is  either  -  Yes  or  No)  to  indicate 
availability  of  that  type  of  data. 

7.2.1  Meta-Rules  at  Level  1. 

Meta-Rules  at  this  level  correspond  to  Fig.  7.1  and  are  concerned  with  requests 
for  activation  of  modules  for:  control  selection,  control  computation 
procedures,  procedures  that  estimate  control  accuracy,  and  resolution  of 
conflicting  or  ambiguous  information. 

Some  of  the  possible  conclusions  at  this  level  are: 

Cl.l:  The  required  control  accuracy  cannot  be  achieved. 

Cl. 2:  The  required  control  accuracy  can  be  achieved. 

The  first  rule  is  a  negative  result.  It  simply  states  that  if  all  availeJale 
control  for  SAR  imagery  has  been  tried  without  success,  then  there  is  nothing 
more  to  be  tried. 

Rule  1.1; 

If  (sensor  -  SAR)  and  (GPS  -  No)  and 
(AM( INS, planimetry, shadows, shading, E)  -No)),  then  (Cl.l) 
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7.2.2  Meta-Rules  at  Level  2. 


Meta-rules  at  this  level  correspond  to  selections  at  the  level  of  Fig.  7.2,  and 
are  concerned  with  selecting  the  most  appropriate  category  of  control  from  the 
available  categories. 

Hie  first  rule  concerns  the  superior  accuracy  of  GPS. 

Rule  2.1: 


If  (GPS  -  Yes),  then  (CT(GPS)) 

Hie  next  rule  accepts  the  high  relative  accuracy  of  INS,  but  requires  some 
additional  control  for  achieving  absolute  accuracy.  A  hierarchy  of  desirable 
control  supplements  to  IMS  is  now  given.  Hie  first  choice  will  be  the  use  of 
pl2unimetric  data. 

However,  the  appropriate  planimetric  data  must  be  available  and  identified 
before  control  computation  and  evaluation  can  take  place.  Hie  precise 
identification  of  particular  planimetric  data  will  take  place  by  the  backward¬ 
chaining  of  rules  in  the  lower  levels  of  the  hierarchy.  But  at  this  level,  all 
that  is  required  is  an  indication  of  whether  or  not  such  data  is  available,  euid 
if  so,  its  usage  in  the  functions  AM  and  CT. 

Rule  2.2: 


If  ((GPS  -  No)  and  (INS  -  Yes)  and  (planimetry  -  Yes),  then 
(CT( INS, planimetry)  and  AM( INS, planimetry, E) ) 

Rule  2.3: 


If  ((H2.2  «  Yes)  and  (AM(INS,  planimetry, E)  -  Yes),  then  (Cl. 2) 

Hie  next  rule  chooses  terrain-induced  shadow  matching  as  a  supplement  for 
planimetry  in  the  case  of  SAR  imagery. 

Rule  2.4: 


If  ((H2.2  -  Yes)  and  (AM( INS, planimetry, E)  -  No)  and  (sensor  -  SAR)  and 
(shadow  -  Yes)),  then  (CTUNS,  planimetry,  shadow)  and 
(AM( INS, planimetry, shadow, E) ) ) 

Rule  2.5: 

If  ((H2.4  -  Yes)  and  (AM( INS, planimetry, shadow,E)  -  Yes),  then  (Cl. 2) 
Hie  next  rule  chooses  terrain-induced  shading  to  supplement  shadow  matching. 

Rule  2.6: 


If  ((H2.3  -  Yes)  and  (AM( INS, planimetry, shadow, E)  -  No))  and  (shading  - 
Yes)),  then  (CT(INS, planimetry, shadow, shading) )  and 
(AM(lNS,planimetry, shadow, shading) ) ) 
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7.2.3  Meta-Rules  at  Level  3; 


The  rules  at  this  level  correspond  to  selections  within  each  of  the  modules  in 
Fig.  7.2  for  GPS,  INS,  SAR  shadows,  shading,  or  planimetry. 

7. 2. 3.1  Planimetry 

The  rules  for  evaluating  planimetry  are  first  partitioned  in  Fig.  7.3.  There, 
the  methods  involve  basic  physical  principles  and  "rules  of  th\]mb”.  A  candidate 
planimetric  feature  is  separately  evaluated  in  terms  of  each  of  these  two 
methods,  and  a  conpssite  rank  evaluation  is  then  made. 

This  con^josite  evaluation  will  also  be  used  to  determine  whether  any  feature's 
score  is  high  enough  to  merit  the  declaration;  planimetry  -  Yes.  This  means  that 
there  is  at  least  one  feature  vhich  is  reason^ly  capable  of  being  detected  in 
the  imagery  under  the  circumstances. 

7. 2. 3. 1.1  Planimetry  Rules  of  Thumb 

The  heuristics,  or  niles  of  thxanb,  relevant  to  the  extraction  of  planimetry 
involve  the  capedaillties  to  both  sense  features  and  extract  these  sensed 
features  from  imagery,  as  shown  in  Fig.  7.5. 

7. 2. 3. 1.1.1  Image  Formation  Rules  of  Thumb 

Some  of  the  heuristics  that  must  be  incorporated  as  rule-based  knowledge  about 
imaging  features  include  a  list  of  commonly  recognizable  objects.  For  SAR,  these 
include; 

o  natural 

-terrain  relief 


-lakes 

-rivers 


-tree  lines 


-shorelines 


o  man-made 

-power  lines 
-railway  lines 
-pipelines 
-seismic  lines 


-roads 
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-bridges 

-dams 

-urban  areas 
-power  plants 
-airfields 
-strip  mines 
-oil  platforms 

Other  knowledge  that  must  be  incorporated  are  heuristics  concerning  conditions 
for  visibility  and  non-visibility  in  SAR  images  for  such  features  as  those 
listed  eibove.  Such  heuristics  include: 

o  Lakes  may  or  not  create  specular  returns  depending  on  the  surface 
roughness. 

o  Lake  boundaries  may  be  difficult  to  distinguish  precisely  because  of 
water  level  cheuiges  cuid  imaged  mud  beunks. 

o  Corner  reflections  make  urban  areas  visible,  but  with  indistinguishable 
boundaries.  Therefore,  such  regions  are  only  suitable  for  very 
approximate  location  estimation. 

o  Thin  objects,  such  as  pipelines  and  power  lines,  are  sometimes  not 
visible  themselves  but  often  lie  in  visible  forest  clearings. 

o  Roads  are  often  only  visible  if  they  are  roughly  in  the  direction  of 
the  flight  line. 

o  Features  that  should  be  visible  coiild  be  obscured  by  other  nearby,  tall 
features  because  of  radar  layover. 

o  Sometimes,  near  range  "striping"  is  observed  that  is  due  to  sidelobe 
effects. 

o  A  false,  linear  signature,  punctuated  by  periodic  highlights,  can  be 
due  to  imaging  a  moving  aircraft. 

It  will  be  important  to  incorporate  these  and  other  "rules  of  thumb"  for 
evaluating  feature  matches.  They  will  also  be  importemt  for  giving  possible 
expleuiations  of  why  certain  features  which  should  ordinarily  be  detectable  in 
computed  image  locations,  based  on  other  feature  matches,  are  not  visible  in  a 
particular  situation. 
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8.0  Spatial  Database  Issues 

The  issues  addressed  in  this  section  concern  appropriate  data  structures  for 
spatial  databases,  and  the  database  management  fmctions  for  interfacing  the 
control  generation  p' -  'tedures  with  a  spatial  database. 

8.1  Data  Structures 

The  following  discussion  represents  a  selective  overview  of  data  structures 
suitable  for  internally  representing  spatial  data  in  a  database.  Such  structures 
must  be  organized  to  efficiently  support  the  geometrically-oriented  queries  that 
are  common  to  spatial  database  transactions.  Typical  examples  of  such 
geometrical  operations  are  containment,  closest  point,  adjacency,  eund 
intersection  [NcKeown,84] .  Additionally,  such  a  database  must  support  the 
maintenance  of  various  source  types  of  spatial  data,  such  as  imagery,  maps, 
terrain  models,  etc. 

A  number  of  such  structures  representing  spatial  data  have  been  developed.  One 
approach  is  via  spatial  decomposition.  Another  approach  encodes  feature 
information  in  the  form  of  point  coordinates  and  vectors. 

Recursive  deconposition  of  the  area  or  volume  to  be  represented  is  the  basis  for 
methods  enploying  quadtrees,  ninetrees,  and  segments.  Hexagonal  structures  have 
been  used  for  non-recursive  decomposition. 

Reviews  of  the  quadtree  concept  can  be  foiand  in  [Samet,84].  Although  quadtrees 
in  their  basis  form  are  strictly  areal  structures,  techniques  have  been 
developed  to  encompass  the  representation  of  linear  features  and  the  boundaries 
of  areal  features.  Such  techniques  involve  conversion  algorithms  to  place 
boundary  data  or  vector  datas  into  a  quadtree  representation  and  vice-versa 
[Samet,84],  [Dyer, 80],  [Nelson, 86].  The  use  of  "linear  quadtrees" 
[Gargantini,82]  increases  the  efficiency  with  v4iich  quadtree  data  can  be 
manipulated  within  memory.  Another  generalization  of  the  quadtree  concept  is  k-d 
trees  [Bently,75],  vrtiich  are  based  on  the  correspondence  between  quadtree 
representations  for  point  data  and  multidimensional  generalizations  of  binary 
search  trees  [Finkel,  74]. 

Recursive  segmenting  decon^sition  of  linear  features  into  smaller  eund  smaller 
segments,  each  one  represented  by  an  elongated  rectangle  containing  it,  is  the 
basis  of  strip  trees  [Ballard, 81] . 

Three-by-three  areal  decomposition  is  the  basis  of  the  nine-tree  concept 
[Raetzsch,85] .  This  structure  offers  some  computational  simplifications  over  the 
quadtree . 

Hexagonal  decomposition  leads  to  an  addressing  scheme  known  as  the  Generalized 
Balauiced  Ternary  representation  [Gibson, 84]. 

Data  structures  based  on  point  and  vector  coordinate  lists  have  been  examined  in 
[Kropatsch,81] .  In  this  formulation,  point  coordinates  are  maintained  in  lists 
ordered  so  that  the  vectors  and  polygons  defined  by  those  points  are  represented 
by  the  location  of  those  coordinates  in  the  list.  Additionally,  pointers  are 
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used  to  link  related  vectors  represented  in  differing  places  in  the  coordinate 
list.  A  geographical  information  system  being  develop^  at  Kork  Systems,  Inc. 
uses  such  a  point/vector  format  with  associated  attributes  in  a  layered 
organization  in  order  to  group  features  by  location  and  size.  The  design  of 
query  languages  vrfiich  efficiently  exploit  such  data  structures  for  efficient 
search  is  considered  in  [ Frank, 1982 ] . 

Ihe  greatest  advantage  of  a  point/\'ector  format  for  map  database  management  auid 
the  generation  of  ground  control  information  for  registration  is  the  logical 
similarity  between  this  type  of  representation  and  the  likely  ground  control 
features  present  in  maps  and  images.  Such  features  are  usually  line  segments, 
curved  boundaries,  intersections  of  segments  and  bounrlaries,  and  centroids  of 
areal  features.  Areal  decon^sition  representations,  such  as  quadtrees  and 
similar  structures,  are  not  as  well  suited  to  the  representation  of  such 
features . 

8.2  Database  Management  Functions 

In  general,  the  management  of  spatial  data  imposes  special  requirements  on  a 
database  management  system  (DBMS).  Moreover,  the  particular  application  of 
control  generation  must  be  considered.  The  type  of  spatial  DBMS  which  is 
relevant  to  this  section  is  one  which  supports  the  particular  application  of 
control  information  generation. 

To  this  end,  spatial  databases  must  support  geometrically  and  spatially  oriented 
queries,  as  well  as  queries  which  involve  non-geometric,  relational  entities. 

Some  typical  queries  to  spatial  databases  that  are  relevant  to  the  image  control 
problem  addressed  in  this  effort  involve  computations  such  as:  geographical 
containment  as  a  function  of  location  xincertainties,  line-of  sight  obscuration, 
existence  and  usage  of  sufficient  terrain  relief,  as  well  as  feature  type, 
location,  and  direction.  Moreover  queries  may  often  occur  as  combinations  of  the 
above. 

A  traditional  relational  organization  could,  in  principle  support  such  non¬ 
relational  database  queries.  However,  non-spatially  oriented  relational 
databases  have  structural  and  performance  problems  associated  with  handling 
spatial  data.  Among  them  are  the  variable  length  records  associated  with 
varicU^le  length  lines,  the  dispersion  of  attributes  over  multiple  teibles,  and  a 
generally  awkward  link  structure  between  features  [ Keating, 87 ] . 

Therefore,  to  achieve  efficient  performaunce,  the  organizational  structure  of  the 
database  should  be  conpatible  with  the  different  types  of  queries  that  are 
expected. 

This  implies  that  a  spatial  database  should  contain  multiple  levels  of 
orgauiization  to  support  the  different  types  of  queries.  Therefore,  data  must  be 
organized  on  a  geographic  and  geometric  basis,  as  well  as  the  traditional 
relational  basis. 


Further  improvements  in  performance  can  be  achieved  if  the  physical 
organization  of  geometric  data  is  closely  compatible  with  its  logical 
organization.  In  particular,  disk  file  records  for  geographically  clustered 
items  should  also  be  physically  clustered.  Additionally,  memory  resident  DBMS 
operations  and  internal  buffering  according  to  entity  types  and  spatial  domains 
will  also  increase  response  time  efficiencies  [Keating, 87] . 


A  factor  v^ich  influences  flexibility  in  the  continued  use  of  such  databases  is 
object-oriented  design.  Such  an  approach  would  more  easily  allow  the  addition 
and  deletion  of  attributes  of  geometric  and  geographic  data  without  requiring 
changes  in  the  overall  logical  structure  of  the  dateibase. 
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9.0  Sumnary  and  Conclusions 


Some  conclusions  can  be  dravm  from  this  preliminary  vrork.  A  terrain  database  ceui 
be  exploited  to  obtain  matching  data  between  a  real  and  synthetic  SAR  image  for 
the  purposes  of  obtaining  image  control.  However,  the  resection  in  space  will 
generally  be  ill-conditioned  if  all  the  matches  are  within  a  narrow  value  for 
ramge.  Therefore,  what  is  required  is  rugged  terrain  features  for  casting  radar 
shadows  at  more  than  one  rauige  offset,  or  the  use  of  another  source  of  control, 
such  as  feature  matching,  at  other  remge  offsets.  These  issues  are  discussed  in 
sec.  3. 

Feature  matching  also  has  good  potential  for  confuting  image  control,  provided 
that  suitable  features  can  be  selected,  emd  that  an  efficient  search  process  is 
used  to  locate  them  them  in  the  image.  These  issues  are  discussed  in  sec.  3  and 
7. 
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Itils  section  briefly  describes  some  of  the  general  Issues  Involved  in 
calibrating  imagery  for  registration  with  maps  or  for  stereo.  Especially 
important  is  the  problem  of  resection  in  space,  or  exterior  orientation. 


10.1  Photoqrammetry 

Most  frame-based  optical  camera  systems  are  based  on  central  perspective 
transformations  of  3-D  Euclidean  space  to  2-D  Euclidean  space.  Other  sensor 
types,  such  as  Synthetic  Aperture  Radar  [Leberl,83]  or  scanned  image  arrays 
require  a  different  geometric  model.  We  will  assume  central  perspective  as  the 
Idealized  transformation  model  in  the  following  discussion. 

The  full  determination  of  a  stereo  model  also  requires  knowledge  of  the  relative 
geometries  of  the  two  stereo  images,  a  model  of  the  distortions  vdiich  are 
perturbations  to  the  central  perspective  model,  and  the  relationship  of  the 
imaging  sensor  to  a  ground  coordinate  system.  In  photogrammetry,  such 
relationships  are  compris^  in  the  determination  of  four  types  of  orientation: 
interior,  exterior,  relative,  and  absolute. 

For  a  single  image,  the  interior  and  exterior  orientations  are  relevant. 
Determining  the  interior  orientation  refers  to  the  process  of  reconstructing  the 
internal  geometry  of  light  bundles  within  the  optical  system.  The  relevant 
parameters  to  be  determined  are  the  focal  length,  the  image  coordinates  of  the 
principal  point,  and  a  model  of  the  lens  distortions  v^ich  cause  deviation  from 
a  central  perspective  transformation.  Typical  of  such  distortions  are  symmetric 
radial  distortion  and  asymmetric  lens  decentering.  Coefficients  quantifying 
these  parameters  are  obtained  during  the  camera  calibration  process. 

Determining  an  exterior  orientation  refers  to  the  process  of  relating  the 
coordinate  system  of  the  imaging  sensor  to  a  grotind  coordinate  system.  Ihis 
essentially  requires  the  determination  of  the  geographic  position  of  the 
exposure  center  and  the  direction  of  the  optical  2ucis. 

For  a  pair  of  images,  the  relevant  orientations,  neglecting  distortions,  are  the 
relative  and  absolute  orientations.  In  principle,  one  could  determine  eui 
interior  and  an  exterior  orientation  for  each  image  separately,  and  from  these 
also  obtain  the  relative  relationship  between  both  images.  However,  such  a 
procedure  would  not  deal  with  measurement  errors  and  inconsistencies  as  well  as 
using  both  images  simultaneously  to  estimate  the  relative  geometry  between  them 
and  to  relate  image  coordinates  to  ground  coordinates.  Distortions  must  be 
accounted  for  by  either  an  interior  orientation  separately  for  both  images,  or 
in  a  bundle  adjustment  for  both  images  simultaneously. 

The  relative  orientation  of  a  stereo  pair  is  concerned  with  the  relative 
position  emd  attitude  of  the  image  coordinate  systems  of  the  two  images.  The  aim 
of  a  relative  orientation  determination  is  to  find  that  relative  orientation  of 
both  images  such  that  all  light  bundles  through  corresponding  image  points 
intersect.  Because  of  the  inevitable  distortions  in  the  lenses,  the  atmosphere, 
and  photographic  processing,  a  "best  fit"  orientation  is  determined  which 
minimizes  some  total  measure  of  bundle  intersection  discrepancy. 


The  parameters  that  are  sought  are  two  relative  position  coordinates 
perpendicular  to  the  line  segment  joining  the  perspective  centers,  and  three 
rotational  transformation  parameters  describing  the  rotation  between  the  two 
coordinate  systems.  The  length  of  the  segment  joining  the  perspective  centers  is 
usiially  not  computed  in  the  determination  of  relative  orientation,  since  this 
parameter  influences  the  scale  of  the  model.  Scale  is  determined  in  absolute 
orientation. 

Absolute  orientation  determination  involves  determining  the  scale  of  the  model, 
"leveling”  the  model,  and  determining  the  position  of  its  origin  in  a  ground 
coordinate  system. 

Bundle  adjustment  is  a  method  for  integrating  the  two  operations  of  relative 
amd  absolute  photo-orientations  into  a  single  operation,  possibly  using  a  number 
of  images  vdiich  cover  some  area.  This  method  can  be  used  for  simultaneously 
estimating  distortions  in  a  stereo  pair  of  images.  Such  a  simultemeous 
estimation  method  allows  for  better  estimates  thaui  are  availe±>le  from  individual 
interior  orientation  determinations. 

In  the  case  of  one  pair  of  images,  this  process  would  decompose  the  problem  into 
a  non-statistical  and  a  statistical  subproblem.  The  non-statistical  portion 
would  involve  determining  an  arbitrary  scale  and  orientation  for  the  model. 
These  can  be  accomplished  by  arbitrarily  assigning  (x,y,z)  coordinates  to  two 
corresponding  points  in  the  image  pair,  emd  by  specifying  a  model  rotation  about 
the  line  connecting  these  two  points.  These  two  points  and  the  rotational  angle 
constitute  seven  degrees  of  freedom. 

The  statistical  portion  would  involve  identifying  some  number  of  corresponding 
point  pairs  to  estimate  the  lens  distortion  coefficients. 

10.2  Radargramimetry 

"Resection  in  space"  is  the  photogrananetric  term  for  calculation  of  the  elements 
of  the  exterior  orientation  (see  section  9.1)  for  a  central  perspective  imaging 
sensor.  This  orientation  is  the  set  of  geometric  pareuneters  that  describe  the 
imaging  process.  For  a  central  perspective  optical  camera,  it  involves  the 
position  of  the  perspective  center  and  the  orientation  of  the  focal  plane. 

For  a  SAR,  it  is  a  description  of  the  flight  path,  during  imaging,  as  a  fxinction 
of  time.  Such  a  flight  path  for  an  aircraft  would  ideally  be  a  straight  line  of 
constant  altitude.  However,  because  of  aerodynamic  effects,  this  is  obviously 
not  achievable  in  practice.  Short-term  departures  from  such  a  straight-line 
flight  path  can  be  estimated  from  the  INS  or  INS/GPS  data.  Especially  important 
are  the  velocity  errors  and  distance  offsets  from  the  nominal  flight  line. 

Using  such  data,  real  time  phase-history  and  range-gating  corrections  can  be 
implemented  in  a  real  time  SAR  processor  to  process  the  data  as  if  it  were  being 
imaged  along  a  nominal  straight  line  flight  path  with  constant  velocity  (see 
section  2.3.4).  Residual  errors  still  exist,  of  course,  because  of  imperfect 
corrections. 
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Because  of  t±iese  residual  errors,  however,  the  radargrainmetric  formulation  of 
resection  in  space  requires  a  time-dependent  formulation  of  the  flight  path 
parameters.  These  range-projection  equations  are  [Leberl,83): 


r  -  Ip  -  ¥(t)|  -  0  (range  equation) 

s(t).(p  -  s(t))  -  sinX  |s(t)|  Ip  -  s(t)|  -  0  (Doppler  equation) 
vdiere: 


p  is  a  control  point 

r  is  rzmge  from  sensor  to  control  point  at  time  t 

lit) ,  s(t)  are  the  instantaneous  position  amd  velocity  vectors  of 
the  sensor  platform 

X  is  the  sensor  squint  angle  (squint  usually  zero  for  aircraft 
SAR,  nonzero  for  satellite  SAR) 

Itiese  nonlinear  equations  are  usually  linearized  using  a  Taylor's  series 
expaunsion.  The  system  of  equations  is  then  solved  for  the  unknown  linear 
correction  for  the  original  estimate.  This  process  is  repeated  iteratively. 
Convergence  requires  good  initial  estimates. 

Af  described  in  (Raggam,87],  problems  of  accuracy  and  convergence  remain  because 
of  discrepancies  between  the  actual  imaging  scenario  and  the  initial  describing 
parameters  for  imaging.  Because  of  these  initial  inaccuracies,  modified 
procedures  have  been  implemented  which  attempt  to  refine  the  estimates  of 
parameters  with  the  greatest  inaccuracies: 

o  range  offset  (constant  for  azimuthal  lines) 

o  azimuthal  shift  (caused  by  time  calibration  errors) 

o  scaling  factors  for  pixel  spacings 

Some  problems  with  initial  estimates  also  arise  because  of  correlations  among 
parameters.  This  seems  to  be  especially  a  problem  in  the  orbital  SAR  case,,  and 
can  involve  [Raggam,87]; 

o  azimuth  shift  and  squint 

o  range  offset,  pixel  spacings,  and  sensor  position 

For  rugged  terrain  imaged  by  satellite  SAR  imagery,  it  was  found  in  [Raggam,87] 
that  siraultemeous  estimation  of  all  parameters  was  not  possible  in  the  orbital 
case,  and  that  step-wise  refinements  were  required. 
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